Identification of Prognostic Biomarkers and Correlation With Immune Infiltrates in Hepatocellular Carcinoma Based on a Competing Endogenous RNA Network

Abstract

Background: Hepatocellular carcinoma (HCC) is one of the most common malignant tumors worldwide. Recently, competing endogenous RNAs (ceRNA) have revealed a significant role in the progression of HCC. Herein, we aimed to construct a ceRNA network to identify potential biomarkers and illustrate its correlation with immune infiltration in HCC.

Methods: RNA sequencing data and clinical traits of HCC patients were downloaded from TCGA. The limma R package was used to identify differentially expressed (DE) RNAs. The predicted prognostic model was established using univariate and multivariate Cox regression. A K-M curve, TISIDB and GEPIA website were utilized for survival analysis. Functional annotation was determined using Enrichr and Reactome. Protein-to-protein network analysis was implemented using SRTNG and Cytoscape. Hub gene expression was validated by quantitative polymerase chain reaction, Oncomine and the Hunan Protein Atlas database. Immune infiltration was analyzed by TIMMER, and Drugbank was exploited to identify bioactive compounds.

Results: The predicted model that was established revealed significant efficacy with 3- and 5-years of the area under ROC at 0.804 and 0.744, respectively. Eleven DEmiRNAs were screened out by a K-M survival analysis. Then, we constructed a ceRNA network, including 56 DElncRNAs, 6 DEmiRNAs, and 28 DEmRNAs. The 28 DEmRNAs were enriched in cancer-related pathways, for example, the TNF signaling pathway. Moreover, six hub genes, CEP55, DEPDC1, KIF23, CLSPN, MYBL2, and RACGAP1, were all overexpressed in HCC tissues and independently correlated with survival rate. Furthermore, expression of hub genes was related to immune cell infiltration in HCC, including B cells, CD8⁺ T cells, CD4⁺ T cells, monocytes, macrophages, neutrophils, and dendritic cells.

Conclusion: The findings from this study demonstrate that CEP55, DEPDC1, KIF23, CLSPN, MYBL2, and RACGAP1 are closely associated with prognosis and immune infiltration, representing potential therapeutic targets or prognostic biomarkers in HCC.

Keywords: biomarkers; competing endogenous RNA network; hepatocellular carcinoma; immune infiltration; prognostic prediction model.

FIGURE 1

The flowchart of the present research. HCC: hepatocellular carcinoma. TCGA: The Cancer Genome Atlas. PPI: protein–protein interaction. GO: Gene Oncology. KEGG: Kyoto Encyclopedia of Genes and Genomes. ceRNA: competing endogenous RNAs. LASSO: the least absolute shrinkage and selection operator.

FIGURE 2

The hierarchical clustering heatmaps and volcano plots for all screened differentially expressed mRNA, miRNA and lncRNA in HCC based on TCGA data. Heatmaps located in the left panels represent differential expressed (DE) mRNAs (A), miRNAs (C), and lncRNAs (E). Volcano plots located in the right panels indicate DEmRNAs (B), DEmiRNAs (D), and DElncRNAs (F) with the cutoff criteria of fold change ≥ 2 and P-value < 0.05. Red color: upregulated, green: downregulated, gray: not statistic expressed.

FIGURE 3

The prognostic predictive model for HCC. (A) The analyzed result of multivariate Cox proportional hazards regression involved with recurrence in TCGA HCC cohort. The middle point of the line indicates the hazard ratio (HR), and the whole length on behalf of the 95% CI for each DEmiRNA. (B) Nomogram based on differentially expressed miRNAs to predict survival in HCC asymptomatic individuals. The prognostic model aim to estimate the survival rate for individual patient, meanwhile reveal the upregulated or downregulated type for each miRNA. At first, draw a line straight upwards from each miRNA to obtain the points from the points axis. Repeat this step until the total scores were gained for 23 miRNAs. Then, after calculating the overall points according to the total points axis, draw a line straight down to the 3-year and 5-year survival axis from the location of total point axis based on the obtained overall scores to indicate the rate for the specific patient (for e.g., the 3-year survival rate is 60% if a patient get the total points of 400).

FIGURE 4

The assessment of miRNAs-based clinical prediction model. The calibration curves according to nomogram model to estimate the survival rate at 3-year (A) and 5-year (B), the X and Y-axis represent predicted and actual survival time respectively. The efficacy of prognostic model of 3- and 5-years survival rate based on time dependent receiver operated characteristic curves (C). The Kaplan-Meier curve of overall survival time between the high- and low- risk groups stratified by the mean of total risk scores (D).

FIGURE 5

The lncRNA-miRNA-mRNA regulatory network in HCC cohort visualized by Cytoscape software 3.6.1 (A). Rectangle represent the 6 DEMs, diamond represent the 56 DELs, the ellipse represent the 28 DEGs. The functional enrichment analysis of DEGs correlated to ceRNA network (B–F). Top 10 biological process (BP) terms (B). Top 10 cell components (CC) terms (C). Top 10 molecular functions (MF) terms (D). Top 10 significantly KEGG pathways (E). Top 15 enriched Reactome pathways (F).

FIGURE 6

The construction of hub gene associated with ceRNA network based on analysis of protein to protein interaction (PPI) network. The PPI network of DEGs (A). PPI network of 6 hub genes, the color of nodes from red to yellow indicates that the connected degrees between each molecule with others decrease gradually (B). the hub gene ceRNA regulatory network including 3 DEMs (rectangle) and 6 hub DEGs (ellipse) as well as 46 lncRNA (diamond) (C).

FIGURE 7

The transcriptional expression of six hub genes in HCC (right panels) and normal liver samples (left panels) based on Oncomine database (A–F). The relative mRNA expression of six hub genes in HCC cell lines and liver cell line (G–L).*P < 0.05, **P < 0.01. The protein expression of hub genes based on Human Protein Atlas database (M–N). The protein expressed level of DEPDC1 was not detected in both tumor and normal liver tissues and the expressed data of CLSPN was lacking in database.

FIGURE 8

(A–F) Six hub genes expression correlated with macrophage polarization in HCC. Markers include CD86 and CSF1R of monocytes; CD68 and IL10 of TAMs (tumor-associated macrophages); IRF5 and PTGS2 of M1 macrophages; and CD163, VSIG4, and MS4A4A of M2 macrophages.