Bioinformatics analysis of the potentially functional circRNA-miRNA-mRNA network in breast cancer

Abstract

Breast cancer (BC) is the most common cancer among women with high morbidity and mortality. Therefore, new research is still needed for biomarker detection. GSE101124 and GSE182471 datasets were obtained from the Gene Expression Omnibus (GEO) database to evaluate differentially expressed circular RNAs (circRNAs). The Cancer Genome Atlas (TCGA) and Molecular Taxonomy of Breast Cancer International Consortium (METABRIC) databases were used to identify the significantly dysregulated microRNAs (miRNAs) and genes considering the Prediction Analysis of Microarray classification (PAM50). The circRNA-miRNA-mRNA relationship was investigated using the Cancer-Specific CircRNA, miRDB, miRTarBase, and miRWalk databases. The circRNA-miRNA-mRNA regulatory network was annotated using Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database. The protein-protein interaction network was constructed by the STRING database and visualized by the Cytoscape tool. Then, raw miRNA data and genes were filtered using some selection criteria according to a specific expression level in PAM50 subgroups. A bottleneck method was utilized to obtain highly interacted hub genes using cytoHubba Cytoscape plugin. The Disease-Free Survival and Overall Survival analysis were performed for these hub genes, which are detected within the miRNA and circRNA axis in our study. We identified three circRNAs, three miRNAs, and eighteen candidate target genes that may play an important role in BC. In addition, it has been determined that these molecules can be useful in the classification of BC, especially in determining the basal-like breast cancer (BLBC) subtype. We conclude that hsa_circ_0000515/miR-486-5p/SDC1 axis may be an important biomarker candidate in distinguishing patients in the BLBC subgroup of BC.

Fig 1. The steps of the BC PAM50 subtype analysis.

The miRNA–mRNA interactions were estimated with mirDB (v6), miRTarBase (Release 8.0), and miRWalk (v3). DE: Differentially expressed, CSCD (v2.0): The Cancer-specific circRNAs database, BC: Breast cancer, FC: Fold change, log: logarithm base 2.

Fig 2. The volcano plot for DECs in BC based on the two microarray datasets from GEO and intersected up- and down regulated circRNAs.

The volcano plot for DEMs in BC based on the two microarray datasets from TCGA and EGA, and intersected up- and down regulated miRNAs. (A): GSE101124, (B): GSE182471, (C): The intersected up- and down-regulated circRNAs between the GSE101124 and the GSE182471 datasets. (D): TCGA, (E): METABRIC, (F): The intersected up- and down-regulated miRNAs from shared miRNAs in the TCGA and the METABRIC datasets, The intersection of the down-regulated mRNAs (G) and the up-regulated mRNAs (H) between PAM50 subtypes, DECs: differently expressed circRNAs, BC: Breast cancer, DEMs: Differentially expressed miRNAs, EGA: European Genome-phenome Archive, BC: Breast cancer, hsa: Homo-sapiens.

Fig 3. The combined violin and box plots for the normalized expression values of hsa_circRNA_101004, hsa_circRNA_000585, hsa_circRNA_100435 in GSE101124 and GSE182471 datasets by the unpaired two-samples Wilcoxon test according to tumor and control samples.

(A): hsa_circRNA_101004 in GSE101124, (B): hsa_circRNA_000585 in GSE101124, (C): hsa_circRNA_100435 in GSE101124, (D): hsa_circRNA_101004 in GSE182471, (E): hsa_circRNA_000585 in GSE182471, (F): hsa_circRNA_100435 in GSE182471.

Fig 4. The combined violin and box plots for the normalized expression values of miR-141-5p, miR-183-5p, and miR-486-5p in TCGA and METABRIC datasets by the unpaired two-samples t-test according to basal and control samples.

(A): miR-141-5p in TCGA, (B): miR-183-5p in TCGA, (C): miR-486-5p in TCGA, (D): miR-141-5p in METABRIC, (E): miR-183-5p in METABRIC, (F): miR-486-5p in METABRIC.

Fig 5. circRNA–miRNA–mRNA regulatory network.

The network consisting of three cricRNAs (hsa_circRNA_000585, hsa_circRNA_101004, and hsa_circRNA_100435), three miRNAs (miR-486-5p, miR-141-5p, and miR-183-5p) and 339 genes was generated by Cytoscape 3.9.0.

Fig 6. Log2FC values of the selected 18 DEGs.

Fig 7

The correlation heatmap of the selected mRNA and miRNAs (A):All PAM50 groups, (B): Only BLBC subgroup, Survival analysis of SDC1. (C): Overall survival of SDC1 in TCGA dataset, (D): Disease-free survival of SDC1 in TCGA dataset, (E): Disease-free survival of SDC1 in GSE25066 dataset.

Fig 8. Identification of hub genes from the PPI network by bottleneck algorithm using the cytoHubba Cytoscape plugin.

The node color changes gradually from blue to red in ascending order according to the log2 (fold change) of genes. (A): A PPI network of the 339 target genes playing crucial roles in BC. This network consists of 138 nodes and 780 edges. The node size changes gradually from small to large in ascending order according to the number of the PMIDs from DisGeNET per gene. (B): A PPI network consist of the ten hub genes (colored blue and red) and 4 extended genes (colored gray) extracted from a. This network consists of 14 (10 hub genes and 4 extended genes) nodes and 24 (14 between hub genes and 10 between extended genes) edges. PPI protein–protein interaction, BC: Breast Cancer.

Fig 9. CircRNA–miRNA–hubgene network.

The network consisting of three circRNAs (hsa_circRNA_000585, hsa_circRNA_101004, and hsa_circRNA_100435), three miRNAs (miR-486-5p, miR-141-5p, and miR-183-5p) and 10 hub genes (AHNAK, CAV1, CDK1, EGR1, FGF2, FOS, KIF11, PPARG, SDC1, and TNXB) was generated by Cytoscape 3.9.0.

Fig 10

Top five Gene Ontology (GO) enrichment annotations of the ten hub genes: (A): biological process, (B): cellular component, (C): molecular function. (D)The significantly enriched Kyoto Encyclopedia of hub-genes and genomes (KEGG) pathways with a FDR < 0.05. The results of the GO and the KEGG analyses were obtained from the ‘Enrichr’ web tool (https://maayanlab.cloud/Enrichr/) and visualized by R package ‘ggplot2’. Cohort plot shows that the ten hub genes are correlated via ribbons with their assigned KEGG terms. FDR: False discovery rate, is calculated using the Benjamini-Hochberg method to adjust the multiple hypothesis testing.

Fig 11. The summary of the possible role of circRNA/miRNA/gene axis in BC pathogenesis according to our study.