Large-scale analyses identify a cluster of novel long noncoding RNAs as potential competitive endogenous RNAs in progression of hepatocellular carcinoma

Differential and clustering analysis

The workflow chart was shown in Supplementary Figure 1. We first identified 11036 differentially expressed genes (DEGs) from three profiles, GSE29721, GSE40367, and GSE62232, consisting of 3949 (35.8%) upregulated and 7087 (64.2%) downregulated genes. The top 5 DEGs in the HCC samples compared with normal samples from different profiles were shown in Table 1 according to their fold change (FC) values. Besides, 3826 lncRNAs were selected as differentially expressed lncRNAs (DELs) from the same three profiles, including 615 (16.1%) upregulated and 3211 (83.9%) downregulated lncRNAs. From another two profiles, GSE36915 and GSE74618, 206 differentially expressed miRNAs (DEMs) were selected, including 80 (38.8%) upregulated and 126 (61.2%) downregulated miRNAs. The top 10 DELs and DEMs in the HCC samples in contrast to normal samples from different profiles were shown in Supplementary Tables 1 and 2, respectively. The results of hierarchical cluster analyses were shown in heatmaps for the expression level changes of DEGs (Figure 1A), DELs (Figure 1B), and DEMs (Figure 1C) between HCC and normal tissues.

Figure 1. Cluster analysis of differentially expressed profiles. Hierarchical cluster dendrogram of DEGs (A) and DELs (B) identified in GES29721, GSE40367 and GSE62232 as well as DEMs (C) identified in GSE36915 and GSE74618. The rows showed DEGs, DELs, and DEMs, while the columns showed paired samples. The pink part represents normal samples and the blue part represents HCC samples. The left vertical axis shows clusters of DEGs, DELs and DEMs, while the above horizontal axis shows clusters of samples. Red represents high expression and green represents low expression.

Intersection analysis

Next, we selected 1016 intersecting DEGs and 116 intersecting DELs among the three profiles to perform further analysis. Twenty-one DEMs were also selected between the two profiles GSE36915 and GSE74618. The intersections were displayed as Venn diagrams (Supplementary Figure 2A–2C).

GO and pathway analysis

According to the intersecting DEGs, we identified 1803 upregulated DEGs-related GO terms and 6339 downregulated DEGs-related GO terms using functional enrichment analysis. The plots of the top 25 upregulated and downregulated GO enrichment terms were shown in Figure 2A and 2B. The results demonstrated that the upregulated DEGs mainly participated in cell division, sister chromatid cohesion, mitotic spindle organization, DNA replication, and mitotic cell cycle; while the downregulated DEGs were closely associated with the oxidation-reduction process, xenobiotic metabolic process, and the epoxygenase P450 pathway.

Figure 2. Top25 enrichment of GO and KEGG pathway analyses for upregulated and downregulated DEGs. The upregulated (A) and downregulated (B) DEGs enriched in GO categories. The upregulated (C) and downregulated (D) DEGs enriched in different pathways. The horizontal axis represents the enrichment score of DEGs. The vertical axis represents different GO categories and pathways. The bubble size indicates the number of genes in each category and pathway, and different colors correspond to different log (FDR) values.

We subsequently distinguished 221 pathways among the upregulated DEGs and 274 pathways among the downregulated DEGs using KEGG pathway analysis based on the intersecting DEGs. The plots of the top 25 pathways for the upregulated and downregulated DEGs were shown in Figure 2C and 2D. The most obviously upregulated pathways were the cell cycle, the p53 signaling pathway, and DNA replication. Pathways dramatically enriched in downregulated genes were metabolic pathways, chemical carcinogenesis, and retinol metabolism. Taken together, these dramatically changed GO terms and pathways might be involved in HCC initiation and progression.

Signal-net analysis

To clarify the interactions among different genes and their products, we further constructed a gene signal-net using 1016 intersecting DEGs based on KEGG database according to the network biology theory. Signal-net could display the relationships between different gene groups and identify upstream and downstream molecules by obtaining gene interactions in multiple pathways. The network we constructed included 97 interacting DEGs, consisting of 23 upregulated genes and 74 downregulated genes (Supplementary Figure 3). According to their regulation degree, several DEGs were considered as the hub genes that exerted the most significant regulatory function in the network. For example, PRKAA2, which held the highest degree, had the strongest interactions with other genes in the network. The pathway information from the KEGG database indicated that PRKAA2 and its related genes were involved in the AMPK/PI3K/AKT signaling pathways, a representative signaling pathway involved in the development of HCC [16]. Besides, several closely connected hub genes, including ITGA6, ITGA2, and ITGB3, are the main cellular adhesion receptors belonging to the integrins family and extensively implicated in multiple steps from cancer initiation to metastasis by acting as signaling molecules, mechanotransducers, and critical components of the cell migration machinery [17]. Therefore, DEGs with a high regulation degree in this signal-net, including PRKAA2, PLCB1, and several genes encoding the integrins family members might play crucial roles in the regulation of HCC development.

CeRNA network

The above results demonstrated that some intersecting DEGs might play a crucial role in HCC development by participating in important processes, such as cell cycle, p53 signaling pathway, and metabolic pathways. We next investigated the underlying molecular mechanisms regulating these processes. LncRNAs have been reported to promote tumor progression through various types of gene regulatory mechanisms, such as epigenetic and transcriptional regulation and serving as ceRNAs for miRNAs [10, 11, 14, 15, 18–20]. Therefore, we constructed a ceRNA regulatory network to reveal the unknown mechanisms driving HCC development.

Intersection datasets were acquired between the DEGs involved in the significant enriched GO terms and KEGG pathways with P < 0.05 and false discovery rate (FDR) < 0.05, which included 393 DEGs. We predicted 213 target genes and 85 target lncRNAs that might be regulated by the 21 intersected DEMs among the 393 DEGs and the 116 intersecting DELs, respectively. The key lncRNA-miRNA and miRNA-mRNA pairs were shown in Supplementary Tables 3 and 4, respectively. According to the constructed association and the theory of ceRNAs, we chose the negatively correlated lncRNA-miRNA and miRNA-mRNA pairs to build an mRNA-miRNA-lncRNA network using Cytoscape v3.0 (Figure 3). As a result, 59 lncRNAs were identified as ceRNAs that interacted with 17 essential miRNAs, and subsequently indirectly regulated 63 coding mRNAs involved in HCC development.

Relying on the constructed relationships in the ceRNA network, we inferred that lncRNAs might indirectly participate in several significant KEGG pathways by serving as ceRNAs of mRNAs, such as the PI3K-AKT signaling pathway enriched by ITGA2, PRKAA2, and CDK1, the cell cycle enriched by MCM5, CDC6, and CHEK1, as well as ECM-receptor interactions enriched by ITGA6 and ITGA2. Similarly, we also inferred that lncRNAs participated in some significant biological processes, including cell division enriched by BIRC5 and CDCA5, as well as cell proliferation enriched by CD34 and ITGA2. Indeed, some miRNA-mRNA and lncRNA-miRNA pairs in this network have been verified to play a part in cancer development in previous studies. For instance, Jafarzadeh et al. provided experimental evidences for hsa-miR-497-5p as a negative regulator of SMAD3, which was a key modulator of the TGF-β signaling pathway during carcinogenesis [21]. Among the genes interacting with hsa-miR-497-5p in this network, BIRC5, a well-known cancer-related gene encoding survivin, has been proven to be upregulated by TGF-β to modulate the cell cycle and apoptosis in various types of cancer [22, 23]. Therefore, the regulatory relationship between hsa-miR-497-5p and BIRC5 presented in the ceRNA network was justifiable. Besides, for lncRNAs, MCM3AP-AS1 has been reported to directly bind to miR-194-5p and act as ceRNA, which subsequently facilitated the expression of miR-194-5p’s target gene FOXA1 in HCC cells, thus promoting HCC cell proliferation, colony formation, and cell cycle progression [24]. Another lncRNA identified in this network, HAND2-AS1, has been proven to increase cell migration of HCC cell lines [25].

Figure 3. The ceRNA network. (A) Network constructed by upregulated miRNAs, downregulated lncRNAs and downregulated mRNAs. (B) Network constructed by downregulated miRNAs, upregulated lncRNAs and upregulated mRNAs. Red diamonds represent upregulated miRNAs, red balls, upregulated mRNAs, and red triangles, upregulated lncRNAs. Green diamonds represent downregulated miRNAs, green balls, downregulated mRNAs, green triangles, downregulated lncRNAs.

PPI (protein-protein interaction) network analysis

Subsequently, we built a protein-protein interaction network predicting the interaction among the proteins encoded by the 51 DEGs in the ceRNA network (Supplementary Figure 4). In the network, ESR1, IGF1, BIRC5, and CD34 had higher degrees (17, 15, 12, and 10, respectively) (Supplementary Table 5). The genes encoding these proteins have been confirmed to be associated with HCC progression [26–29]. Combined with the functional analysis results, we found that the interacting DEGs in the ceRNA network were mainly enriched in cell proliferation, cell adhesion, the PI3K-Akt signaling pathway, the p53 signaling pathway, pathways in cancer, and metabolic pathways.

Survival analysis and expression validation for the ceRNA network in TCGA dataset

In the Gene Expression Omnibus (GEO) datasets, we identified an extensive and comprehensive lncRNA-miRNA-mRNA regulation network in HCC development. To further clarify their expression and prognostic value, we performed Kaplan-Meier survival analysis for all the DEGs, DELs, and DEMs in the ceRNA network in patients with HCC from TCGA. The results showed that 6 lncRNAs (Figure 4), 23 mRNAs, and 4 miRNAs had a significant impact on OS (Supplementary Figure 5) and a consistent expression pattern (Supplementary Figure 6) in the GEO datasets. Among them, 3 lncRNAs (GBAP1, MCM3AP-AS1, and SLC16A1-AS1) (Figure 4A–4C, 4G–4I), 15 mRNAs (ACLY, BAX, BIRC5, CDC6, CDCA5, CHEK1, FOXM1, ITGA2, LAMC1, MCM5, MYB, PAK1, PRKAA2, SOX4, and TK1) (Supplementary Figures 5Aa–5Ao, 6Aa–6Ao) and 3 miRNAs (hsa-miR-10b-5p, hsa-miR-183-5p, and hsa-miR-222-3p) (Supplementary Figures 5Aa–Ao, 6Ba–Bc) were identified as pro-tumor factors because of their high expression in cancer tissues and their correlation with shorter OS in patients with HCC. In contrast, another three lncRNAs (C3P1, DIO3OS, and HNF4A-AS1) (Figure 4D–4F, 4J–4L), eight mRNAs (ABAT, ACSM2A, ASPA, CAMK4, CYP8B1, ESR1, IGF1, and PDE7B) (Supplementary Figure 5Ca–5Ch, 6Ca–6Ch) and one miRNA (hsa-let-7c-5p) (Supplementary Figures 5D, 6D) showed low expression in cancer tissues and correlated with longer OS, implying that these transcripts might be protective factors in HCC.

Figure 4. The expression and survival significance of the six key lncRNAs in the ceRNA network in HCC TCGA database. Kaplan-Meier survival curves showed significant OS differences between high- and low-expression of GBAP1 (A), MCM3AP-AS1 (B), SLC16A1-AS1 (C), C3P1(D), HNF4A-AS1 (E) and DIO3OS (F) in TCGA dtabase. Expression levels of GBAP1 (G), MCM3AP-AS1 (H), SLC16A1-AS1 (I), C3P1(D), HNF4A-AS1 (E) and DIO3OS (F) in normal and HCC tissues from TCGA database. N, normal tissue; Ca, cancer tissue.

According to the interactions in the ceRNA network, we inferred that the aberrant expression of DEGs might be indirectly regulated by six lncRNAs (GBAP1, MCM3AP-AS1, SLC16A1-AS1, C3P1, DIO3OS, and HNF4A-AS1). To further determine the association of the identified lncRNAs with survival time, we conducted univariate and multivariate Cox regression model analyses based on the clinical characteristics in TCGA dataset. As shown in Table 2, in the univariate analysis, all six lncRNAs could be incorporated in the COX regression model. Multivariate analysis demonstrated that the expression levels of lncRNAs C3P1, DIO3OS, and SLC16A1-AS1 were independent prognostic factors for OS in patients with HCC. Thus, it was reasonable to infer that these lncRNAs might be crucial factors to predict the prognosis of patients with HCC.

Correlation analysis in TCGA dataset

Since the main aim of this study was identifying the clinical noteworthy lncRNAs and uncovering the regulatory role of lncRNAs as potential ceRNAs to mediate downstream/upstream RNAs functions/ associations, we selected the six lncRNAs with great survival significance in the ceRNA network for further study. To validate the regulatory role of the six lncRNAs as ceRNAs, we analyzed the correlation with their associated miRNAs and mRNAs in the ceRNA network based on TCGA data containing 366 patients with HCC. As shown in Figure 5, a considerable number of lncRNA-miRNA pairs showed a significant negative correlation for both upregulated (Figure 5A) and downregulated (Figure 5D) lncRNAs. Some of the six lncRNAs-associated miRNA-mRNA pairs also demonstrated a negative correlation (Figure 5B and 5E). Moreover, almost all six lncRNAs-associated lncRNA-mRNA pairs, for either upregulated (Figure 5C) or downregulated (Figure 5F) lncRNAs, demonstrated a tight positive correlation in TCGA database. Together, these results indicated the proposed six lncRNAs, GBAP1, MCM3AP-AS1, SLC16A1-AS1 C3P1, DIO3OS, and HNF4A-AS1, had a partially negative correlation with their associated miRNAs, as well as a partially positive correlation with their associated mRNAs, further confirming that the six might function as ceRNAs to regulate the expression of mRNAs and miRNAs in the progression of HCC.

Figure 5. Pearson correlation analysis for six key lncRNAs with their associated miRNA and mRNAs in 366 HCC patients from TCGA. Pearson correlograms of lncRNA-miRNA paris (A), lncRNA-mRNA pairs (B) and mRNA-miRNA pairs (C) in upregulated lncRNAs including GBAP1, MCM3AP-AS1 and SLC16A1-AS1. lncRNA-miRNA paris (D), lncRNA-mRNA paris (E) and mRNA-miRNA paris (F) in downregulated lncRNAs including C3P1, HNF4A-AS1 and DIO3OS by Pearson correlation analysis. The correlation coefficient R value ranges from -1 to 1, the color of which changes from blue to brown.

Expression validation and survival analysis for the six lncRNAs in clinical samples

To further verify the significance of the six lncRNAs mentioned in section 2.8, we analyzed their expression levels in tumor tissues and adjacent non-tumor tissues from 158 diagnosed patients with HCC using quantitative real-time reverse transcriptase-polymerase chain reaction (RT-PCR). As shown in Figure 6A–6F, all six lncRNAs were differentially expressed in tumor tissues and normal tissues. Three lncRNAs (GBAP1, MCM3AP-AS1, and SLC16A1-AS1) were upregulated in HCC tissues (Figure 6A–6C), while the other three lncRNAs (C3P1, HNF4A-AS1, and DIO3OS) were downregulated in HCC tissues (Figure 6D–6F). These results were consistent with the previous bioinformatic analysis. Furthermore, based on the clinical information, we analyzed the relationship between the expression levels of these DELs and OS, progression-free survival (PFS) and distant metastasis-free survival (DmFS) in these patients, which was demonstrated using Kaplan-Meier curves in Figure 6G–6X. Consistent with the results in TCGA, the upregulated lncRNAs, including GBAP1, MCM3AP-AS1, and SLC16A1-AS1, were associated with worse OS (Figure 6G–6I), PFS (Figure 6M–6O) and DmFS (S-U), thus identifying them as pro-tumor factors. Inversely, the downregulated lncRNAs, including C3P1, HNF4A-AS1, and DIO3OS, were correlated with the better OS (Figure 6J–6L), PFS (Figure 6P–6R) and DmFS (Figure 6V–6X), and thus represented protective factors. Importantly, combined with the clinical characteristics, we constructed univariate and multivariate Cox proportional hazards regression models for OS, PFS, and DmFS. As shown in Table 3, multivariate Cox’s regression analysis revealed that DIO3OS, MCM3AP-AS1, and SLC16A1-AS1 expression levels were independent prognostic factors for OS. C3P1, MCM3AP-AS1, and SLC16A1-AS1 expression levels were independent prognostic factors for PFS, and HNF4A-AS1, MCM3AP-AS1, GBAP1, and SLC16A1-AS1 expression levels were independent prognostic factors for DmFS. These data indicated that the proposed six lncRNAs were clinical noteworthy biomarkers for predicting the prognosis and metastasis for patients with HCC.

Figure 6. Expression and survival sifnificance for the six lncRNAs in 158 clinical samples. Expression levels of lncRNAs GBAP1 (A), MCM3AP-AS1 (B), SLC16A1-AS1 (C), C3P1 (D), HNF4A-AS1 (E), and DIO3OS (F) in 158 HCC tissues and adjacent normal tissues. Kaplan-Meier survival curves showed significant OS differences between high- and low-expression of lncRNAs including GBAP1 (G), MCM3AP-AS1 (H), SLC16A1-AS1 (I), C3P1(J), HNF4A-AS1 (K), and DIO3OS (L), and significant PFS differences including GBAP1 (M), MCM3AP-AS1 (N), SLC16A1-AS1 (O), C3P1 (P), HNF4A-AS1 (Q) and DIO3OS (R), as well as significant DmFS differences including GBAP1 (S), MCM3AP-AS1 (T), SLC16A1-AS1 (U), C3P1 (V), HNF4A-AS1 (W) and DIO3OS (X). N =158, Ca = 158; N, normal tissue; Ca, cancer tissue.

Correlation validation for the six lncRNAs in clinical samples

To further understand the links of the six lncRNAs with miRNAs or mRNAs in the ceRNA network, we selected lncRNAs associated miRNAs, as well as specific mRNAs with survival significance in TCGA dataset of HCC from the ceRNA network. And their expression levels were assessed by RT-PCR in 158 clinical samples. By correlation analysis based on RT-PCR results, a significantly negative correlation was found between a majority of lncRNA-miRNA pairs for both upregulated (Figure 7A) and downregulated (Figure 7D) lncRNAs. In accordance with TCGA dataset, most of the lncRNAs-associated miRNA-mRNA pairs demonstrated a negative correlation in clinical samples (Figure 7B and 7E), and the positive correlation in a great part of lncRNA-mRNA pairs was also observed (Figure 7C and 7F). Furthermore, we analyzed the expression pattern of the proposed ceRNA network in samples from different stages (early = TNM I or II; late = III) based on the results of RT-PCR, including the six lncRNAs and associated miRNAs and specific mRNAs. As shown in Figure 7G, we found the upregulated lncRNA GBAP1, MCM3AP-AS1 and SLC16A1-AS1 were highly expressed in late-stage samples compared with early-stage samples. Inverse and similar expression patterns were observed in their associated miRNAs and mRNAs, respectively (Figure 7G). In contrast, the downregulated lncRNA C3P1, DIO3OS, and HNF4A-AS1 were low expressed in late-stage samples compared to early-stage samples (Figure 7H). Inverse and similar expression patterns were also observed in their associated miRNAs and mRNAs, respectively (Figure 7H). Taken together, these data further verified that the six lncRNAs might play an important role as ceRNAs in regulating the expression of miRNAs and mRNAs. Moreover, upregulation of GBAP1, MCM3AP-AS1, and SLC16A1-AS1 and downregulation of C3P1, DIO3OS, and HNF4A-AS1 might be involved in regulating the progression of HCC from early to late stage.

Figure 7. Correlation validation for the six lncRNAs with their associated miRNA and mRNAs in 158 clinical samples of HCC. RT-PCR was performed to detect the expression of the six lncRNAs associated miRNA and mRNAs in 158 clinical samples. Pearson correlation analysis for six key lncRNAs with their associated miRNA and mRNAs was performed based on RT-PCR results. Pearson correlograms of lncRNA-miRNA paris (A), lncRNA-mRNA pairs (B) and mRNA-miRNA pairs (C) in upregulated lncRNAs including GBAP1, MCM3AP-AS1 and SLC16A1-AS1. lncRNA-miRNA paris (D), lncRNA-mRNA paris (E) and mRNA-miRNA paris (F) in downregulated lncRNAs including C3P1, HNF4A-AS1 and DIO3OS by Pearson correlation analysis. The correlation coefficient R value ranges from -1 to 1, the color of which changes from blue to yellow. 158 clinical samples were diveded into early and late groups according to TMN stage (early = TNM I or II; late = III), the expression of the six lncRNAs associated miRNA and mRNAs was shown as heatmaps. (G) Expression of the upregulated lncRNA GBAP1, MCM3AP-AS1 and SLC16A1-AS1 and their ssociated miRNA and mRNAs. (H) Expression of th edownregulated lncRNA C3P1, DIO3OS, and HNF4A-AS1 and their ssociated miRNA and mRNAs. Log₂(P) value was shown as color ranging from blue to red.

Receiver operating characteristic (ROC) curves analysis for the key lncRNAs in clinical samples

To assess the discriminatory ability of the six key lncRNAs in HCC, ROC curve analyses were conducted in 158 patients with HCC and the areas under the curve (AUCs) were calculated. As shown in Figure 8A–8F, in the assessment of OS, the AUCs of five lncRNAs (GBAP1, MCM3AP-AS1, SLC16A1-AS1, C3P1, and DIO3OS) were more than 0.75, and the AUCs of HNF4A-AS1 was 0.737. In the assessment of PFS, the AUCs of all the six lncRNAs were greater than 0.75 (Figure 8G–8L). Additionally, five lncRNAs (GBAP1, MCM3AP-AS1, SLC16A1-AS1, HNF4A-AS1 and DIO3OS) had a good performance in diagnosing DmFS with AUCs above 0.75, and the AUCs of C3P1 was 0.712 (Figure 8M–8R). These results suggested that the six lncRNAs had good sensitivity and specificity to predict survival and distant metastasis in patients with HCC.

Figure 8. ROC curves for the six key lncRNAs in 158 clinical samples. ROC curves of lncRNA GBAP1, MCM3AP-AS1, SLC16A1-AS1, C3P1, HNF4A-AS1 and DIO3OS for OS (A–F), PFS (G–L), and DmFS (M–R) respectively. The AUCs under binomial exact confidence interval was calculated to generate the ROC curve.

Patients and samples

Tissue samples, including HCC tumor tissues and adjacent non-cancerous tissues (n = 158), were obtained from patients who had resection of primary HCC in the Cancer Center of Sun Yat-sen University between 2005 and 2008. None of these patients had received preoperative chemotherapy or radiotherapy. After resection, matched fresh tissues were immersed immediately in RNAlater® (Ambion, Austin, TX, USA), kept overnight at 4 °C, and then stored at −80 °C until RNA isolation before qRT-PCR detection. Follow-up was performed in our outpatient department and involved clinical and laboratory examinations every three months for the first two years, every six months during the third to fifth years, and annually for an additional five years or until death, whichever occurred first. Follow-up periods for survivors ranged from 2 to 87 months, with a median follow-up of 41 months. OS and PFS were used as measures of prognosis. Written informed consent was obtained from each patient, and the Ethics Committee of Sun Yat-Sen University Cancer Center approved the study protocol.

GEO gene expression datasets

HCC gene expression data were obtained from the GEO database (http://www.ncbi.nlm.nih.gov/geo/), including GSE29721, GSE40367, and GSE62232. The three profiles included a total of 120 samples, consisting of 98 cancerous and 22 normal samples. The mRNA expression profiles were obtained using an Affymetrix Human Genome U133 Plus 2.0 Array. According to the annotation file, we chose the probe sets and BLAST searched them in the NCBI RefSeq database to identify the probe sets of noncoding RNAs with a length of over 200 bp. We then used these probe sets to acquire the expression data of lncRNAs from the three profiles. The miRNA expression data were also obtained from the GEO databases. Two profiles, GSE36915 and GSE74618, were selected that contained 286 cancerous and 31 normal samples.

Differential and clustering analysis

Based on the expression data from the datasets, a random variance model (RVM) t-test was applied to filter the DEGs, DELs, and DEMs among cancerous and normal samples. After significance analysis and FDR analysis, we selected the DEGs, and DELs from GSE29721, GSE40367, and GSE62232 according to the P-value threshold. P < 0.05 and FC > 2 was considered as significant difference. DEMs were selected using a threshold of P-value < 0.05, FDR < 0.05, and FC > 1.5. Hierarchical cluster analysis was performed and a cluster dendrogram was constructed to demonstrate distinct characterizations of screened DEGs, DELs, and DEMs between the cancerous and normal tissues. Furthermore, we selected the intersecting DEGs, DELs, and DEMs among the different profiles.

GO and KEGG pathway analysis

Gene functional enrichment analysis was used to predict the biological functions of the intersected DEGs according to the GO database with P < 0.05 and FDR <0.05. The potentially involved signaling pathways were identified using the KEGG pathways analysis program (http://www.genome.jp/kegg/tool/map_pathway1.html)).

Two-side Fisher’s exact test and a χ² test were used to classify the GO categories, and the FDR [47] was calculated to correct the P-value: The smaller the FDR, the smaller the error in judging the P-value. The FDR was defined as:

$$\text{FDR}=\text{1}-N_k/T,$$(1)

where N_k referred to the number of Fisher’s test P-values less than the test P-values. We computed P-values for the GO categories of all the DEGs. Within the significant categories, the enrichment, Re, was calculated as:

$$\text{Re}=\left(n_f\text{/}n\right)\text{/}(N_f\text{/}N),$$(2)

where “n_f” was the number of flagged genes within the particular category, “n” was the total number of genes within the same category, “N_f” was the number of flagged genes in the entire microarray, and “N” was the total number of genes in the microarray [48]. For KEGG pathway analysis, we again used Fisher’s exact test and χ² test to select the significant pathways, and the threshold of significance was defined P-value <0.05 and FDR < 0.05. The Re value was calculated using equation (2) [49- 51].

Signal-net analysis of intersecting DEGs

Based on the KEGG database analysis, a gene-gene interaction network of DEGs was constructed to demonstrate the regulatory relationships among the DEGs identified in the intersection analysis among different expression profiles of mRNAs. The networks were stored and presented as graphs, where nodes were mainly genes (or proteins or compounds) and edges represented the type relationships between the nodes, e.g. activation or phosphorylation. We investigated the nature of networks using tools implemented in the R software.

The important nodes were identified computationally. To this end, we used the connectivity (also known as degree), which was defined as the sum of connection strengths with the other network genes:

$${\text{K}}_i={\displaystyle \sum _{u\ne i}{a}_{ui}}$$,(3)

In the gene networks, the connectivity measured how a gene correlated with all other network genes. For a gene in the network, the number of source genes of a gene was called the indegree of the gene and the number of target genes of a gene was its outdegree. The character of a gene was described using betweenness centrality measures reflecting the importance of a node in a graph relative to other nodes. For a graph G: (V, E) with n vertices, the relative betweenness centrality ${C^{\prime }}_B(v)$ is defined by:

$${C^{\prime }}_B(v)=\frac{2}{n^2-3n+2}{\displaystyle \sum _{\begin{array}{l}s\ne v\ne t\in V\\ s\ne t\end{array}}\frac{{\sigma }_{st}(v)}{{\sigma }_{st}}},$$(4)

where ${\sigma }_{st}$ is the number of shortest paths from s to t, and ${\sigma }_{st}(v)$ is the number of shortest paths from s to t that pass through a vertex v [52-56].

Target transcripts of DEMs prediction and ceRNA network construction

Based on the functional DEGs identified in the GO and KEGG pathway analysis, we constructed the intersection datasets between the DEGs involved in the significant enriched GO terms and pathways with P < 0.05 and FDR < 0.05. Ultimately, 393 DEGs were selected, comprising 78 upregulated DEGs and 315 downregulated DEGs. Combining the intersecting DELs and DEMs, we predicted the targeted sponge lncRNAs using the miRanda tool (http://www.microrna.org/ microrna/home.do), and the target mRNAs using Targetscan (http://www.targetscan.org/) and miRWalk (http://129.206.7.150/). For each pair of miRNA-mRNA or mirRNA-lncRNA, we conducted Pearson correlation analysis and chose the significantly correlated pairs [57]. In particular, we summarized the intersecting targeted mRNAs discerned using miRanda and miRWalk.

Then, we chose the miRNAs that negatively regulated the expression levels of lncRNAs and mRNAs to construct a ceRNA network according to the normalized signal intensity of the expression of specific mRNAs and lncRNAs.

Protein regulation network analysis

The STRING online database tool (https://string-db.org/cgi/input.pl) was used to construct a PPI network of the DEGs identified in the ceRNA network. The interacting pairs with a confidence score greater than 0.4 were considered as significant and were retained. The proteins encoded by the hub genes were screened according to the degree of the nodes. The degree represents the number of interaction partners and was calculated using Perl code.

Survival analysis

We screened the HCC RNA-Seq TCGA datasets containing survival information and selected 370 patients with HCC as a dataset to analyze the relationship between the expression level of lncRNAs and mRNAs and OS. Another dataset that contained 371 patients with HCC was selected to analyze the relationship between the expression level of miRNAs and OS. In the clinical samples, we also investigated the correlation of lncRNAs with OS and PFS for 158 patients with HCC. According to the expression level of mRNAs, lncRNAs, and miRNAs in TCGA or our clinical samples, we classified them into two groups: High expression and low expression. Survival curves were displayed using Kaplan–Meier plots. The Wilcoxon log-rank test was used to analyze the survival difference between the high and low groups. All survival analyses in TCGA were conducted using the R package, Survival.

Correlation analysis

To verify the correlation of expression among the miRNAs, lncRNAs, and mRNAs identified in the ceRNA network, we chose the lncRNAs that had significance for the prognosis of patients with HCC and their associated-mRNAs and miRNAs. Pearson correlation analyses were performed among them based on the RNA-Seq data of 366 patients with HCC searched in TCGA database (https://cancergenome.nih.gov/). The correlogram was constructed using the R package, Corrplot.

RNA extraction and qRT-PCR

Total RNA was isolated from tissues using the TRIzol reagent (Invitrogen Corporation, Waltham, MA, USA) according to the manufacturer’s instructions. The concentration and purity of the RNA were evaluated using a NanoDrop 2000 instrument (Thermo Scientific, Waltham, MA, USA). For mRNA, the first-strand cDNA was synthesized from total RNA using a GoScript Reverse Transcription System (Promega, Madison, WI, USA). GAPDH was used as an endogenous control for normalization. For miRNAs, the first-strand cDNA was synthesized from total RNA using a riboSCRIPT^TM Reverse Transcription Kit (RIBOBIO, GuangzZhou, China). U6 was used as an endogenous control for normalization. QPCR was performed using GoTaq qPCR Master Mix (Promega).

ROC curve analysis

ROC analyses were performed using the pROC package in the R language, based on data from 158 clinical samples. The diagnostic ability of the prediction model was evaluated by calculating the area under a ROC curve. The ROC curve was used for classifier evaluation and was drawn by plotting sensitivity against the false-positive rate. The AUC under a binomial exact confidence interval was calculated to generate the ROC curve.

Statistical analysis

Data are shown as the mean ± s.d. and analyzed using Student’s t-test. A paired t-test was used for paired samples. Statistical analyses were performed using GraphPad Prism 7 and SPSS software (GraphPad Software, La Jolla, CA, USA). A value of P<0.05 was considered statistically significant.