Construction of a ceRNA network to reveal a vascular invasion associated prognostic model in hepatocellular carcinoma

2.1 Data acquisition and processing

We downloaded the circRNA expression dataset GSE94508 [17] and GSE97332[18] containing human HCC and adjacent normal tissues from GEO (https://www.ncbi.nlm.nih.gov/geo/) database. Two sets of data are from platform GPL19978 and our search terms were HCC and circRNA. There are 10 samples in GSE94508, which are divided into 5 tumor tissues and 5 tumor adjacent normal tissues, and 14 samples in GSE97332, consisting of 7 tumor tissues and 7 tumor adjacent normal tissues. The primary data were processed by background correction and quantile normalization, and the integrated circRNA expression data were obtained after batch effect removal using R package sva [19]. The ComBat function in R packet sva [19] was used to integrate multiple data and remove the batch effect. We used TCGAbiolinks package [20] to download HCC (Liver hepatocellular carcinoma, TCGA-LIHC)-related gene expression data from the TCGA database (https://portal.gdc.cancer.gov/), which contained 371 HCC tissues and 50 tumor adjacent normal tissues. Meanwhile, TCGAbiolinks package [20] was used to obtain the corresponding clinicopathological survival information of 371 tumor samples, including age, survival status, follow-up time, VI, stage, and so on. In addition, the gene expression data of normal liver tissue (GTEX-Liver) were downloaded from GTEX database, which contained 110 normal samples. Two datasets of TCGA-LIHC and GTEX-Liver were combined to extract the expression data of common genes, and the batch effect between different datasets was corrected by ComBat function in R package sva [19]. The combined data included 371 HCC tumor samples and 160 control samples. Gene expression data included miRNA and mRNA expression levels of each patient, and clinical information included age, sex, pathological stage, VI, survival status, and total survival time. Only patients with complete survival information and gene expression data were included in this study. Finally, a total of 358 patients were selected for further analysis.

2.2 Screening of DEGs

In order to assess the difference of gene expression between HCC and normal tissue, R packet limma [21] was used to analyze the difference among different groups. |logFC| > 1 and P value < 0.05 were set as thresholds for differentially expressed circRNAs (DEcircRNAs), where circRNAs with logFC > 1 were considered significantly up-regulated, and circRNAs with logFC < −1 were considered down-regulated. For differentially expressed mRNAs (DEmRNAs), we set the threshold value as |logFC| > 2 and P value < 0.05, where DEmRNAs with logFC < −2 were down-regulated and DEmRNAs with logFC > 2 were up-regulated. Similarly, |logFC| > 1 and P value < 0.05 were used to identify differentially expressed miRNAs (DEmiRNAs). The volcano map was used to show the up and down-regulated DEGs, and pheatmap R-package [22] was used to draw the heat map of these DEGs in all samples.

2.3 Construction of a ceRNA network

To better understand the effect of circRNAs on the progression of HCC, a ceRNA network was built based on DEcircRNAs, DEmiRNAs, and DEmRNAs. The human sequences of DEcircRNAs and DEmiRNAs were downloaded from circBase (http://www.circbase.org/) and miRBase (version 21; http://www.mirbase.org/) databases, respectively. The miRanda prediction tool was used to predict the interaction between DEcircRNAs and DEmiRNAs. In addition, the potential mRNAs targeted by DEmiRNAs were obtained from StarBase (http://starbase.sysu.edu.cn/) database, which provided prediction results from seven prediction programs (TargetScan, microT, miRmap, picTar, RNA22, PITA, and miRanda). If the interaction between miRNA and mRNA was predicted in not less than four programs, it was selected for next analysis. Then, we overlapped the target mRNAs with the DEmRNAs mentioned above. The nodes that could not form circRNA-miRNA-mRNA interaction relationship were removed, and finally a ceRNA network was established and visualized by the Cytoscape software (version 3.7.0/www.cytoscape.org).

2.4 Functional enrichment analysis

Gene Ontology (GO) functional annotation analysis is a common method for large-scale functional enrichment of genes [23], including biological processes (BPs), molecular functions (MFs), and cellular components (CC). Kyoto Encyclopedia of Genes and Genomes (KEGG) is a widely used database for storing information about genomes, biological pathways, and drugs [24]. In order to study the BPs that the ceRNA network might participate in, we chose the DEmRNAs within the ceRNA network for GO functional annotation analysis and KEGG pathway enrichment analysis by using clusterProfiler R software package [25]. P value < 0.05 was considered statistically significant.

2.5 Gene expression difference between HCC with and without VI

According to the clinical information of patients with HCC, we divided patients into two groups according to HCC with or without VI; 236 patients with HCC marked “None” without VI, and 122 patients with HCC labeled “Micro” and “Macro” with VI. In order to explore the difference between HCC with and without VI, limma R-packet was used to analyze the DEGs between the two groups. |logFC| > 1 and Padj < 0.05 were set as the threshold of DEGs. The genes of logFC > 1 and Padj < 0.05 were regarded as the up-regulated DEGs, and the genes of logFC < −1 and Padj < 0.05 were down-regulated. The volcano map was used to show the up-regulated and down-regulated DEGs, and the pheatmap R-package drew the heat map of these DEGs in all samples.

2.6 Construction of PPI network

PPI network is composed of individual proteins through the interaction between each other, which participates in various life processes such as biological signal transmission, gene expression regulation, and cell cycle regulation. PPI analysis is very important for elucidating the molecular mechanism of key cellular activities in carcinogenesis. Based on the DEGs between HCC with and without VI, the STRING database (https://string-db.org/) was used to evaluate the PPI information and further construct the PPI network [20]. Taking the interaction score of 0.4 as the cut-off value, the PPI network was visualized by cytoscape software.

2.7 Construction of prognostic model and analysis of prognosis

Based on the intersection of DEGs within the PPI and ceRNA network mentioned above, we constructed a prognostic risk score model through the method of least absolute shrinkage and selection operator (LASSO) to realize the risk stratification of HCC patients. LASSO regression is a commonly used variable selection method when fitting high-dimensional generalized linear models. The glmnet package [26] in R was used to execute the LASSO algorithm. The risk scoring model was established by combining the regression coefficient with the corresponding gene expression value. After calculating the risk score, taking the median risk score as the cut-off point, the patients in each cohort were divided into low- and high-risk group accordingly. Univariate and multivariate Cox analyses were used to analyze the ability of risk score in combination with clinicopathological features to predict the OS. Then, the risk score model and clinicopathological parameters were used for further analysis through the clinical predictive Nomogram, which was constructed by rms R-packet [27]. Decision curve analysis (DCA) was performed, and clinical impact curves were drawn to evaluate whether the model-based decisions were beneficial to patients.

2.8 Difference of immune cell infiltration between high- and low-risk groups

The immune microenvironment is mainly composed of immune cells, inflammatory cells, fibroblasts, endothelial cells, bone marrow-derived cells, extracellular matrix, cytokines, and chemokines, which is a comprehensive system. The analysis of immune cells infiltration in the clinical samples is of great importance to discern the mechanisms underlying cancer progression and predict prognosis. Single sample gene cluster enrichment analysis (ssGSEA) is an extension of GSEA method. ssGSEA algorithm was used to calculate the content of 28 kinds of immune cells in high- and low-risk groups [28]. The composition of immune cells in high and low risk groups was visualized by box map.

2.9 Statistical analysis

All statistical analyses were carried out in R language (https://www.r-project.org version 4.0.2). For the comparison of continuous variables in two groups, the statistical significance of normal distribution variables was estimated by independent t-test, the differences between non-normal distribution variables were analyzed by Wilcoxon rank sum test, and the differences between multiple groups of independent variables were analyzed by Kruskal−Wallis Test. The Kaplan−Meier (KM) method was used to evaluate the difference in survival time of patients with HCC, and the logarithmic rank test was used to determine the statistical significance of the observed differences between distinct groups. The hazard ratio and 95% confidence interval were calculated based on Cox regression analysis. All the statistical P values were bilateral, P < 0.05 was considered statistically significant.

1. Ethics approval and consent to participate: TCGA and GEO both are public databases. The patients involved in the database have given their approval in the original studies. Users can download relevant data for research and publish relevant articles. Our study is based on these open-source data, so there are no ethical issues.

3.1 Identification of DEGs

The schematic diagram of our analysis strategy is shown in Figure A1. To analyze the difference of gene expression between HCC and normal tissue, we first analyzed the mRNA expression data. We used limma differential analysis to get 6,674 DEmRNAs, including 4,446 up-regulated DEmRNAs and 2,228 down-regulated DEmRNAs. Using DEmRNAs and data grouping information to draw a classification heat map, DEmRNAs could well distinguish HCC from normal tissue (Figure 1a and b). For miRNA expression data, we used limma differential analysis to get 100 DEmiRNAs, including 33 up-regulated DEmiRNAs and 67 down-regulated DEmiRNAs. Using DEmiRNAs and data grouping information to draw a classification heat map, DEmiRNAs could well differentiate HCC from normal tissue (Figure 1c and d). For the integrated circRNA expression data, we used limma difference analysis to get 211 up-regulated DEcircRNAs and 47 down-regulated DEcircRNAs. Using DEcircRNAs and data grouping information to draw a classification heat map, DEcircRNAs could also well discriminate HCC from normal tissue (Figure 1e and f).

Figure 1

Screening of DEGs between HCC and normal tissue. (a) and (b) Identification of DEmRNA. (c) and (d) Identification of DEmiRNAs. (e) and (f) Identification of DEcircRNAs. Abscissa was log2 FoldChange, ordinate was log10 (adjust P-value), red nodes represent up-regulated DEGs, blue nodes represent down-regulated DEGs, and grey nodes represent non-DEGs.

3.2 Construction of the ceRNA network based on HCC-related DEGs

We used the StarBase database to identify the DEmRNAs targeted by DEmiRNAs. A total of 1,059 miRNA-mRNA interactions were predicted, including 51 DEmiRNAs and 316 DEmRNAs (Figure 2a). Using the miRanda prediction tool, 36 circRNA-miRNA interactions were predicted, based on 19 DEcircRNAs and 25 DEmiRNAs (Figure 2b). After integrating the circRNA-miRNA interaction with the miRNA-mRNA interaction and removing the nodes that cannot form circRNA-miRNA-mRNA interaction, a novel HCC-related ceRNA network was established. The network consisted of 11 DEcircRNAs, 12 DEmiRNAs, and 161 DEmRNAs (Figure 2c).

Figure 2

Construction of the ceRNA network. (a) miRNA-mRNA interaction network, green node represents miRNA and blue node represents mRNA. (b) circRNA-miRNA interaction network, green node represents miRNA and blue node represents circRNA. (c) circRNA-miRNA-mRNA interaction network, green node represents miRNA, orange node represents circRNA, and blue node represents mRNA.

3.3 Functional enrichment analysis of DEmRNAs

In order to explore the biological functions of these HCC-related DEmRNAs, we performed the functional enrichment analysis (Figure 3a and Tables 1 and 2) HCC-related DEmRNAs were mainly enriched in BPs associated with axonogenesis, skeletal system development, mesenchyme development, epithelial cell proliferation, and protein kinase B signaling (Figure 3b). At the same time, it was enriched in CCs such as filopodium, transcription regulator complex, neuron to neuron synapse, neuronal cell body, and postsynaptic density (Figure 3c) and MFs such as DNA-binding transcription activator activity, RNA polymerase II-specificity, DNA-binding transcription activator activity, beta-catenin binding, chemorepellent activity, and growth factor binding (Figure 3d). Then, the pathway enrichment analysis was carried out, and the results showed that HCC-related DEmRNAs were enriched in biological pathways such as breast cancer, melanoma, transcriptional mis-regulation in cancer, p53 signaling pathway, cell cycle, glioma, and so on (Figure 3e). Figure 3f shows the most significant enrichment pathway: hsa05224: breast cancer (Figure 3F). Here we found that it was intriguing to observe the enrichment of a pathway for a different cancer type. Therefore, we compared the DEGs between breast cancer and HCC from TCGA database, and found that there were commonalities between these two cancer types (Figure A2).

Figure 3

GO and KEGG enrichment analysis. (a) GO function enrichment analysis, ordinate was log10 (P value), abscissa was GO terms, node color represented z-score. (b)−(d) The top five items of BP, CC, and MF showed that the node size indicated the number of genes contained in the current GO term, the color of the line represents different GO terms, and the node color represents the FC of the gene. (e) KEGG pathway enrichment analysis showed that the node color represents the FC of the gene, and the line color represents different KEGG pathways. (f) Significantly enriched KEGG pathway, hsa05224: Breast cancer.

3.4 Screening hub genes in relation to VI

A total of 748 DEGs were obtained by differential expression analysis between HCC with VI and those without VI, including 625 up-regulated and 123 down-regulated DEGs (Figure 4a). We screened the proteins encoded by these DEGs to construct a PPI network visualized in Cytoscape (Figure 4b). There are 122 DEGs and 278 PPI pairs in the PPI network. The top five genes that interact with other DEGs are FOS and NRXN1 (interacting with 18 DEGs), CCNA2, ESR1, and NTRK2 (interacting with 17 DEGs). We overlapped the DEGs within the PPI and ceRNA network mentioned above to get seven genes as hub genes (Figure 4c). We used the R packet “GOSemSim” to calculate the GO semantic similarity of these seven hub genes. The results showed that there was a high correlation among HOXC10, HOXC11, HOXC8, and HOXD10 (Figure 4d). Furthermore, the gene expression levels of seven hub genes were screened from two datasets of TCGA-LIHC and GTEX-Liver, and visualized by R-packet ggplot2 box map. The results showed that there were significant differences in the expression of seven hub genes between HCC tumor samples and control samples (Figure A3).

Figure 4

PPI network of DEGs between HCC with and without VI. (a) The DEGs between HCC with and without VI were analyzed, the abscissa was log2 fold change, the ordinate was log10 (adjust P-value), the red node was up-regulated DEGs, the blue node was down-regulated DEGs, and the grey node was non-significant DEGs. (b) The PPI network constructed by the encoded proteins of DEGs, and the orange color represents the key genes. (c) The Venn diagram of genes in PPI and DEmRNA in ceRNA network, the orange color represents the DEmRNA in ceRNA network, and green color represents the genes in PPI network. (d) The GO semantic similarity of key genes, the horizontal axis represents similarity, and the vertical axis was the genes.

3.5 Construction of prognostic risk score model based on hub genes

We used R packet “glmnet” to obtain the regression coefficients of these seven hub genes based on LASSO method. Combined with the gene expression levels of these hub genes, we constructed a risk score model to predict the prognosis of patients with HCC. The regression coefficients of the seven hub genes were as follows: HOXC8: 0.037, GNAO1: 0.017, ADRA1D: 0.020, ARPP21: −0.008, HOXC11: 0.020, HOXC10: −0.014, and HOXD10: 0.047 (Figure 5a and b). We explored the correlation among the expression of these hub genes, and the results showed high levels of correlation among HOXC10, HOXC11, and HOXC8 (Figure 5c). We calculated the risk score for each patient in the TCGA-HCC cohort, and then divided the patients into high and low risk groups according to the median risk score. KM survival analysis showed that OS in the higher risk group was significantly shorter than that in the lower risk group (P < 0.0001, Figure 5d). Next we assessed the relationship between risk score and HCC patient survival. The results showed that the higher the risk score, the shorter the survival time of HCC patient and the higher the proportion of death (Figure 5e).

Figure 5

Construction of risk scoring model. (a) LASSO regression characteristic fitting curve. (b) Quantitative analysis of key genes in LASSO regression. (c) Correlation analysis of expression levels of hub genes. (d) Survival analysis between patients in high- and low-risk groups. Orange color indicated patients in low-risk group, green color indicated patients in high-risk group, the horizontal axis was survival time, and the vertical axis was survival probability. (e) Correlation diagram of the risk factors in risk scoring model. The vertical axis in the top figure was the risk score, the horizontal axis was the patient, the vertical axis in the middle figure was the survival time of the patient, and the horizontal axis was the patient. The lower figure is the heat map of the expression levels of key genes, orange color indicats high expression and green color indicates low expression.

In order to verify the robustness of the risk score, we tested the utility of our model in the validation set data. The verification set data were from the Chinese HCC patients with HBV infection (CHCC-HBV) cohort from 2010 to 2014 at Zhongshan Hospital [29], which was accessed through NODE (https://www.biosino.org/node). Combined with the survival data of patients and our model, the risk score of each patient was calculated, and the patients were grouped and analyzed by risk score. The results showed that the patients with high-risk score in the verification set also showed a trend of worse prognosis (Figure A4).

3.6 Correlation between risk score and clinical parameters

We analyzed the correlation between risk score and clinical parameters, and the results showed that patients in the high-risk group were usually associated with higher T/N/M stages and VI occurrence (Figure 6a). Then, the effects of the expression levels of seven key genes on the prognosis of patients was analyzed, and the results indicated that the high expression of HOXC8, ADRA1D, HOXC11, HOXC10, and HOXD10 predicted a poor prognosis (Figure 6b−f). In order to assess the effects of risk score and clinical characteristics on the prognosis of HCC patients, the Nomogram model was constructed using clinical parameters and risk score to predict the 1, 3, and 5 year survival rate of patients with HCC. Comparing with the actual 1, 3, and 5 year survival rate of patients, the results showed that the survival rate of HCC patients could be well predicted by the Nomogram model (Figure 7b−d). Univariate Cox analysis using risk score and clinical characteristics showed that risk stratification, T/N/M stage, and ethnicity were all prognostic factors for HCC (Figure 7e). Multivariate Cox analysis showed that risk stratification was the only prognostic factor for HCC (Figure 7f). The DCA curve showed that the model could be beneficial to the prognosis evaluation of patients (Figure 7g). The time-ROC analysis also showed that the 1 year survival prediction performance of the prognostic model was 58%, and both the 3 and 5 year survival prediction performance were 61% (Figure 7h).

Figure 6

The prognostic relevance of seven characteristic gene. (a) The correlation between risk score, clinical parameter, and key gene expression level. (b)−(f) The survival curves show the effects of HOXC8, ADRA1D, HOXC11, HOXC10, and HOXD10 gene expression on the prognosis of patients.

Figure 7

The Nomogram models. (a) Nomograph of HCC patients’ clinical parameters and risk score. (b)–(d) Nomogram model was used to predict the 1 year survival rates, 3 year survival rates, and 5 year survival rates of HCC patients. The horizontal axis and vertical axis were the predicted survival rate and the actual survival rate of HCC patients, respectively. (e) and (f) Univariate and multivariate COX analyses of risk score and HCC patients’ clinical parameters. (g) DCA curve was applied to assess the useful degree of the model for HCC patients. (h) Time ROC analysis was applied to predict 2 year survival rates, 3 year survival rates, and 5 year survival rates of HCC patients.

3.7 Differences in immune infiltration

After analyzing the immune cell infiltration in HCC, we found that the counts of activated CD4+ T cell, T follicular helper cell, mast cell, type 2T helper cell, activated CD8+ T cell, regulatory T cell, and natural killer T cell in high-risk group were significantly more than that of low-risk group (Figure 8a). Moreover, the correlations between the immune cells content in the high-risk group and the low-risk group were further calculated, and the results showed that the correlation between the immune cells content in the high-risk group was significantly higher than that in the low-risk group (Figure 8b and c). The relationship between hub gene expression and the immune cells of two groups were analyzed and the result showed that there was a noteworthy difference between HOXC11, HOXC8 expression, and immune cells (Figure 8d−e).

Figure 8

Immune infiltration analysis. (a) The difference of immune cell infiltration between high-risk group and low-risk group, green color indicates high-risk group, orange color indicates low-risk group. (b) and (c) There was a correlation between the immune cell infiltration in the high-risk group and the low-risk group, with blue color indicating a positive correlation and red indicating a negative correlation. (d) and (e) The correlation analysis between immune cell infiltration and hub genes in high-risk group and low-risk group shows that the horizontal axis is immune cell, the vertical axis is key gene, and the node size indicated significant level. The node color represented the correlation level ratio.