Eighteen HCC-related DUcircRNAs were identified
In the GEO database, three sets of HCC-related data were obtained: GSE78520 (three sets of HCC and normal liver tissue), GSE94508 (five sets of HCC and matched paracancerous tissues), and GSE97332 (seven sets of HCC and matched adjacent normal tissues). Analysis identified 175 DUcircRNAs in GSE7850, 201 DUcircRNAs in GSE94508, and 453 DUcircRNAs in GSE97332. Comparison of the identified DUcircRNAs revealed twenty DUcircRNAs that were included in all three datasets (Figure 1A, Table 1). Of these DUcircRNAs, only the mechanism of circZFR and hsa-circ-0000673 were previously reported. DUcircRNAs that were previously unreported are hsa-circRNA-103948, hsa-circRNA-104640, hsa-circRNA-101408, hsa-circRNA-102587, hsa-circRNA-102034, hsa-circRNA-101777, hsa-circRNA-001416, hsa-circRNA-102559, hsa-circRNA-102451, hsa-circRNA-104760, hsa-circRNA-101094, hsa-circRNA-101287, hsa-circRNA-101555, hsa-circRNA-400071, hsa-circRNA-102954, hsa-circRNA-100053, hsa-circRNA-104268, and hsa-circRNA-000996. In Table 1, the maternal genes, their position, and their names in circBase are listed. The differential expression of the 18 DUcircRNAs compared to normal tissues is shown by heat map (Figure 1B).
Final DDmiRNA and DUcircRNA-DDmiRNA network
The miRNAs targeted by the set of 18 final DUcircRNAs were predicted by CircInteractome, respectively, and compared with the 80 DDmiRNAs identified in TCGA. Ten miRNAs were found in both analyses (hsa-miR-136-5p, hsa-miR-145-5p, hsa-miR-187-3p, hsa-miR-370-3p, hsa-miR-377-3p, hsa-miR-383-5p, hsa-miR-758-3p, hsa-miR-874-3p, hsa-miR-326, and hsa-miR-375), and these miRNAs are considered the final DDmiRNAs (Figure 2). According to ANOVA analysis of the patient clinical data, hsa-miR-136-5p was associated with sex (P=0.0321) and pathological stage (P=0.0366); hsa-miR-187-3p was associated with sex (P=0.0253); hsa-miR-383-5p was associated with pathological T status (P=0.0340); hsa-miR-326 was related to pathological T stage (P=0.0097); and hsa-miR-375 was related to sex (P=0.00004). Under multivariate analysis, hsa-miR-326 was still closely related to pathological stage (P=0.0047) and pathological T status (P=0.004), and hsa-miR-187-3p was still related with sex (P=0.019) (Table 2).
Among the 18 DUcircRNAs, hsa-circRNA-102587, hsa-circRNA-101777, hsa-circRNA-102559, hsa-circRNA-400071, and hsa-circRNA-102954 targeted miRNAs were not included in the DDmiRNA analysis, and the remaining 13 DUcircRNAs and their targeted 10 DDmiRNAs formed 25 circRNA-miRNA network pairs, representing the DUcircRNA-DDmiRNA network (Figure 3).
Final DUmRNA and ceRNA networks
According to miRDIP, the 10 DDmiRNAs have the potential to target 7310 mRNAs, of which 3115 mRNAs had integrated scores > 0.5 (the highest integrated scores is 1). There were 316 DUmRNAs in the TCGA LIHC data and 169 mRNAs out of the 3115 mRNAs were identified as DU mRNA. These two sets of mRNAs were combined as the final DUmRNAs, and 169 network pairs were formed by combining these final DU mRNAs with the 10 DDmiRNAs, forming a final ceRNA network (circRNA-miRNA-mRNA network) with the DUcircRNA-DDmiRNA network (Figure 4).
Final DUmRNA functional analysis, PPI network, and hub genes
The list of 169 DUmRNAs were loaded to Metascape, and input as species: H. sapiens, and analysis as species: H. sapiens. The main enriched GO terms were: regulation of neuron differentiation, TP53 Regulates Transcription of Genes Involved in G1 Cell Cycle Arrest, retrograde Trans-synaptic signaling by resonance, synaptic signaling, fat cell differentiation, blood vessel morphogenesis, G0 and Early G1, urogenital system development, connective tissue development, integration of energy metabolism, cell morphogenesis involved in differentiation, learning or memory, negative regulation of DNA binding, response to nutrients, positive regulation of cell development, resolution of D-loop structures through Synthesis-Dependent Strand Annealing (SDSA), brain development, and carbohydrate-derivative biosynthetic process. KEGG analysis showed two enriched terms: Small cell lung cancer and GABAergic synapse (Figure 5).
We next used STRING to establish a PPI network, using the Hide Disconnected Nodes in The Network function to fully display the 47 interactions. As shown in Figure 6, in the PPI network, the line color indicates the type of interaction evidence. Default active interaction sources: Textmining, Experiments, Databases, Co-expression, Neighborhood, Gene Fusion, and Co-occurrence. Based on line color and evidence strength, E2F1, H2AFX, TOP2A, and RAD51 are hub genes associated with the HCC DUcircRNA-DDmiRNA-DU mRNA network. GO and KEGG analysis showed that the hub genes of the PPI network mainly participated in the processes of TP53 Regulates Transcription of Genes Involved in G1 Cell Cycle Arrest and Pathways in cancer and blood vessel morphogenesis (Figure 6).
Pathological Stage and Survival analysis of hub genes
Difference analysis was performed using R software on E2F1, H2AFX, TOP2A, and RAD51 in the TCGA LIHC dataset including 371 HCC patient tissue samples and 50 patients with adjacent normal tissues. Boxplot was used by selecting the Match TCGA normal data option in GEPIA (Figure 7). Compared with normal tissue, E2F1 (logFC = 3.61, P = 2.95E-44), H2AFX (logFC = 1.54, P = 9.35E-16), TOP2A (logFC = 3.82, P = 4.56E-44), and RAD51 (logFC= 2.10, P = 2.69E-25) showed high expression. The log2(TPM + 1) values were used as the Y-axis for a pathological stage plot, which showed that E2F1 (F value = 6.43, Pr (>F)=0.0003), H2AFX (F value = 8.89, Pr (>F)=1.09E -05), TOP2A (F value = 8.15, Pr (>F) = 2.96E-05), and RAD51 ((F value = 7.70, Pr (>F) = 5.42E-05) are closely related to HCC stage (Figure 7). Finally, the prognostic significance of these four hub genes in 364 patients was analyzed using GEPIA. As shown in Figure 7, the four hub genes E2F1 (log-rank P = 0.0025, HR = 1.7), H2AFX (log-rank P = 0.0050, HR = 1.6), TOP2A (log-rank P = 0.0028, HR = 1.7), and RAD51 (log-rank P = 0.0054, HR = 1.6) were significantly associated with overall survival (OS) in HCC patients.
DUcircRNA-DDmiRNA-DUmRNA network related to hub genes
The hub gene-associated DUcircRNA-DDmiRNA-DUmRNA network was selected from the final ceRNA network. The network can be divided into two parts, the hsa-circ-101408/hsa-circ-101094/hsa-circ-000996-hsa-miR-136-5p-E2F1 axis, the hsa-circ-100053/hsa-circ101555/hsa-circ-101094/hsa-circ-101408-hsa-miR-145-5p-H2AFX axis, the hsa-circ-101408-hsa-miR-383-5P-TOP2A axis, and the hsa-circ-104760-hsa-miR-758-3p-RAD51 axis (Figure 8).