Identification of Differentially Expressed circRNAs in BC
Raw data from GSE101123 were downloaded from GEO database and then performed differential gene expression analysis by limma package in R. Using P-adjust < 0.05 and |log2FC|≥ 1 as a threshold, 7 differential expressed circRNAs were screened out, of which, 4 were upregulated and 3 were downregulated. The differential expressed circRNAs were displayed in the volcano plot (Fig.1a), and also exhibited in heat map (Fig.1b). The essential characteristics and basic structural of the seven circRNAs were displayed in Table 1 and Fig 2.
Prediction of circRNA-bound targeted miRNA
Since many circRNAs included the binding site of miRNAs, it could exert sponge activity by binding to a specific miRNA, and therefore could suppress miRNA function. We first predicted circRNA-bound miRNAs using CSCD miRNA target online prediction software, and total 441 target miRNAs of the seven circRNAs were predicted. Furthermore, we downloaded raw sequencing data from the official TCGA website. Using FDR<0.05 and |log2FC|≥ 1 as the threshold, we identified a total of 251 differentially expressed miRNAs, among them 64 were downregulated and 187 were upregulated (Fig.3a). We further used the Venn diagram to intersect the target miRNAs and DEmiRNAs, and screened 5 target miRNAs (miR-934, miR-7760, miR-760, miR-4443, miR-4784) that play important roles in BC (Fig.3b). We then analyzed the correlation between the expression of five miRNAs and the differential expressed circRNAs. Finally, four circRNA-miRNA interactions including 3 circRNAs (circ_0028899, circ_0000375, circ_0000376) and 4 miRNAs (miR-934, miR-7760, miR-4443, miR-4784) were determined.
Prediction of miRNA-Targeted mRNAs
We mapped the abovementioned 4 DEmiRNAs into the miRwalk target gene prediction software and searched for the target gene. Total 7168 target genes of the four miRNAs were predicted. In addition, 553 downregulated differentially expressed genes in BC were obtained from TCGA database (FDR<0.05, |log2FC|≥ 2.5) (Fig.4a). We then used the Venn diagram to intersect the target genes and the downregulated differentially expressed genes, and screened 149 target genes (Fig.4b).
CircRNA–miRNA-mRNA network construction
In order to better comprehend the interaction of circRNA, miRNA and mRNA associated with BC, we established a circRNA-miRNA-mRNA network. As shown in Fig.5, this network was composed of 3 circRNA nodes, 4 miRNA nodes, 149 mRNA nodes and 462 edges. This network showed the possible ceRNA relationships among the three differential expressed circRNAs (circ_0028899, circ_0000375, circ_0000376), the four miRNAs (miR-934, miR-7760, miR-4443, miR-4784), and the 149 differentially expressed genes.
GO and KEGG enrichment analysis of the overlapped genes
The 149 overlapped genes were annotated, and these genes were selected for GO functional enrichment analysis. When considering Biological Process (BP), the overlapped genes were primarily enriched in response to acid chemical, regulation of membrane potential, and regulation of lipid metabolic process. In terms of cellular components (CC), the top three enriched items were plasma membrane protein complex, receptor complex, membrane raft. With regards to molecular function (MF), cation transmembrane transporter activity, peptide binding, and amide binding were enriched in the top three places (Fig.6a). As for the pathways, the overlapped genes were enriched in several KEGG pathways, for instance PPAR signaling pathway, PI3K−Akt signaling pathway, and AMPK signaling pathway (Fig.6b).
Formation of PPI network and hub gene selection of hub genes
We constructed the PPI network by using the string database and visualized it with Cytoscape software to explore the interactions among the 149 overlapped genes. As shown in Fig.7a, a total of 91 nodes and 159 edges were filtered following removal of the unconnected nodes in this network. Furthermore, we screened the hub genes from the PPI regulatory network by applying the CytoHubba plugin in cytoscape. The top 6 hub genes (DGAT2, ACSL1, ADIPOQ, LPL, LEP, PCK1) were screened by using MCC method, and one subnetwork was constructed with 6 nodes and 14 edges (Fig.7b).
CircRNA–miRNA-hub genes network construction
Then, we established a circRNA-miRNA-hub gene network to show the link between circRNAs, miRNA and hub genes. As shown in Fig.8, we discovered nine circRNA-miRNA-hub genes regulatory networks, including circ_0028899/miR-934/LPL, circ_0000375/miR-7706/LPL, circ_0000375/miR-7706/DGAT2, circ_0000375/miR-7706/ADIPOQ, circ_0000376/miR-4443/PCK1, circ_0000376/miR-4443/ADIPOQ, circ_0000376/miR-4443/LEP, circ_0000376/miR-4784/ADIPOQ, circ_0000376/miR-4443/ACSL1.
Assessing the overall survival in patients with BC
We further used Kaplan–Meier plotter to perform prognostic analysis of the six hub genes, the survival analysis results showed that low expression of ADIPOQ, LPL, LEP were associated with poor overall survival in patients with BC (Fig. 9), while the other had no statistically significant effect on overall survival in patients with BC.
Validation of Differentially Expressed circRNAs in BC tissue
To verify the reliability of bioinformatics analysis, the expression of the three circRNAs were detected in 30 pairs breast cancer tissues and adjacent tissues. The results indicated that circ_0028899, circ_0000375, and circ_0000376 were significantly down-regulated in breast cancer tissues (P<0.05, Fig.10), which was consistent with those of bioinformatics analysis.