Identification of breast cancer-related circRNAs by analysis of microarray and RNA-sequencing data: An observational study.

BACKGROUND: An increasing number of studies indicate that circular RNAs (circRNAs) participate in tumorigenesis. The aim of this study was to elucidate the regulatory mechanisms of circRNAs in breast cancer based on the construction of the circRNA-related ceRNA network.
METHODS: The expression profiles of circRNAs, miRNAs, and mRNAs were obtained from the Gene Expression Omnibus (GEO) and The Cancer Genome Atlas (TCGA) databases. A ceRNA network was constructed by Cytoscape. The interactions among proteins were analyzed using the STRING database, and hub genes were extracted using the cytoHubba application. The functions of the differentially expressed mRNAs (DEmRNAs) were analyzed by the Kyoto Encyclopedia of Gene and Genomes (KEGG) and the Gene Ontology (GO) database.
RESULTS: In total, 7 differentially expressed circRNAs (DEcircRNAs), 27 differentially expressed miRNAs (DEmiRNAs), and 102 DEmRNAs were selected for the construction of the ceRNA network of breast cancer. We established a protein-protein interaction network and identified 6 hub genes. Then, a circRNA-miRNA-hub gene regulatory module was established based on 2 DEcircRNAs, 2 DEmiRNAs, and 2 DEmRNAs. GO and KEGG pathway analyses indicated the possible association of DEmRNAs with breast cancer onset and progression.
CONCLUSIONS: The circRNA hsa_circ_0000519 is likely critical in the pathogenesis of breast cancer and may serve as a future therapeutic biomarker.

Introduction

Breast cancer is an example of notorious malignancies among women and is one of the most common causes of cancer-related death.[ According to the latest cancer statistics, approximately 2.1 million women worldwide were diagnosed with breast cancer, and approximately 627,000 individuals died from the condition in 2018.[ The treatment approach and prognosis of breast cancer are based on tumor characteristics, patient factors, and response to treatment. Due to the high intratumor heterogeneity of breast cancer, clinical outcomes vary widely, and the 5-year overall survival remains unsatisfactory.[ Therefore, elucidation of the molecular mechanisms of breast cancer and identification of effective potential molecular biomarkers are essential for clinical decision-making.

As novel non-coding RNAs, circular RNAs (circRNAs) regulate eukaryotic gene expression.[ They are formed by back-splicing covalently joined 3′ and 5′ ends, which is distinct from the canonical splicing of linear RNAs.[ Recent reports indicate that circRNAs have important functions in carcinogenesis and showed diagnostic and prognostic value.[ Salmena et al[ proposed the competitive endogenous RNA (ceRNA) hypothesis in 2011, which posits that there is a complicated regulatory network at the post-transcriptional level and that all types of RNA transcripts serve as the natural sponges of miRNAs to restrain their function through sharing at least 1 miRNA response element (MRE). As a ceRNA, circRNA may interact with miRNA to regulate target gene expression and exert a crucial role in tumor initiation as well as progression.

Previous studies have confirmed that the circRNA-miRNA-mRNA regulatory network has a vital role in the pathogenesis and progression of breast cancer.[ For example, overexpression of circTADA2As inhibits the expression of miR-203a-3p and increases the expression of the miR-124 target gene SOCS3, which inhibits cancer cell proliferation, migration, invasion, and clonogenicity.[ Similarly, circAGFG1 was shown to interact with cyclin E1 (CCNE1) through the sponging of miR-195-5p in triple-negative breast cancer.[ However, the roles of circRNAs and the overall regulatory network in contributing to breast cancer development and progression remain unclear.

To better search for the complex ceRNA network associated with breast cancer, bioinformatics methods were used to analyze differentially expressed genes (DEGs) in breast cancer. The ceRNA network was constructed by integrating data from the Gene Expression Omnibus (GEO) and The Cancer Genome Atlas (TCGA) databases. To better understand the underlying mechanisms contributing to pathogenesis, we performed a series of analyses including functional enrichment analyses and protein–protein interaction analyses. These results can provide further insight into the roles of circRNAs in carcinogenesis and highlight new treatment targets or biomarkers for breast cancer.

Materials and methods

Data collection and DEG screening

The circRNA expression profiles (GSE101123) were obtained from the GEO (https://www.ncbi.nlm.nih.gov/geo/) database, including data from 8 breast cancer tissues and 3 normal tissues. The mRNA (1066 breast cancer tissues and 112 normal tissues) and miRNA (1057 breast cancer tissues and 103 normal tissues) expression profiles were obtained from the TCGA database (https://portal.gdc.cancer.gov). Approval by the Ethics Committee was not necessary because all data were collected from publicly available databases (GEO and TCGA). The raw data were processed by background correction and normalization by using the affy package of R/Bioconductor. The Limma and edgeR packages were used to identify DEGs between normal samples and tumor samples. The threshold was set at | log2 fold change (FC) | > 1.0 and adjusted P-value <.05.

Constructing the ceRNA network

On the basis of differently expressed circRNAs (DEcircRNAs) identified from the microarray dataset, the Cancer-Specific CircRNA database (http://gb.whu.edu.cn/CSCD/) was used to screen the target miRNAs of circRNAs. To maximize the data reliability, these target miRNAs were further screened by the differentially expressed miRNAs (DEmiRNAs) obtained from the TCGA database. Then, miRDB, TargetScan, as well as miRTarBase databases were employed to predict specific to miRNAs.[ Only mRNAs recognized by the 3 databases were considered to be candidate targets and were intersected with the identified differently expressed mRNAs (DEmRNAs) to screen out the DEmRNAs targeted by the DEmiRNAs. Finally, according to the predicted relationship of DEcircRNA-DEmiRNA and DEmiRNA-DEmRNA, a co-expression network (circRNA-miRNA-mRNA) was constructed, and visualized by Cytoscape software (version 3.7.0).

Construction of the PPI network

The STRING database provides protein–protein interaction (PPI) based on a variety of information, such as coexpression, co-occurrence, textmining, and fusion.[ To assess the interactions between DEmRNAs, the protein interactions between DEmRNAs were extracted from the STRING V10 database and displayed with cytoscape 3.7.0. We used the cytoHubba application to explore the hub genes of the obtained PPI network.[

Functional enrichment of DEGs

To provide insight into the underlying biological processes and pathways related to the DEGs in the ceRNA network, the Gene Ontology (GO) as well as Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were carried out through the “clusterProfiler” package in R/Bioconductor.[ The threshold for enrichment significance was P-value <.05.

RNA expression and prognosis

When analyzing the RNAs selected from the ceRNA network, breast cancer patients were assigned to different groups (high or low expression) based on the median gene expression value. We utilized Kaplan–Meier Plotter, which is an online database that includes 5143 breast, 2437 lung, 1816 ovarian, 1065 gastric, and 364 liver cancer patients with overall survival (OS), progression-free survival (PFS), and recurrence-free survival (RFS) data, to assess the relationship between the transcription level of RNAs and prognosis among patients with breast carcinoma. Meanwhile, the log rank P-value was calculated. The hazard ratio (HR) with 95% confidence intervals (CI) was also estimated.[

Results

Identification of DEGs

A total of 7 DEcircRNAs, 4 upregulated and 3 downregulated, were identified in the GSE101123 dataset (Fig. 1). The basic characteristics of the 7 circRNAs are listed in Table 1. A total of 236 DEmiRNAs, 160 upregulated and 76 downregulated, and 5376 DEmRNAs, 4198 upregulated and 1178 downregulated, were identified in the TCGA database (Fig. 2A and B).

Figure 1

Heatmap of the 7 differentially expressed circRNAs from the GSE101123 dataset; the x axis represents the samples, and the y axis denotes the differentially expressed circRNAs. The green color represents the downregulated genes, while the red color represents the upregulated genes.

Table 1

Basic characteristics of the 7 differently expressed circRNAs.

Figure 2

Volcano plot of differentially expressed RNAs. Ascending normalized expression level is colored from green to red. A, miRNA from TCGA; B, mRNA from TCGA. TCGA = The Cancer Genome Atlas.

Construction of the ceRNA network

To further examine the underlying mechanism of circRNAs in mediating mRNA based on miRNA, a circRNA-miRNA-mRNA network was established and subsequently observed by Cytoscape v3.6.0. We retrieved data relating to the top 7 DEcircRNAs identified from the microarray dataset from the Cancer-Specific CircRNA Database (CSCD) online database and identified 448 pairs of interacting circRNAs and miRNAs. After intersecting with the DEmiRNAs, only 48 circRNA-miRNA pairs remained, including 7 DEcircRNAs and 34 DEmiRNAs. We further searched for mRNAs targeted by these 34 DEmiRNAs from the miRDB, TargetScan, as well as miRTarBase databases and selected those overlapping with the identified DEmRNAs. Ultimately, a total of 102 DEmRNAs were involved in the ceRNA network, along with 7 circRNAs and 27 miRNAs (Fig. 3).

Figure 3

The ceRNA network of circRNA-miRNA-mRNA in breast cancer. Vs stands for circRNAs, while diamond is indicative of the miRNAs, rounded rectangles represent the mRNAs, and the grey edges represent the circRNA-miRNA-mRNA interactions. The red nodes indicate upregulated expression, whereas the green nodes indicate downregulated expression.

After removing unconnected nodes, the PPI network was conducted, involving 49 nodes and 55 edges (Fig. 4A). To explore the hub genes in the network, indicating a critical role in the process of breast cancer carcinogenesis, the closeness centrality of DEmRNAs was calculated using the cytoHubba plugin, and the top 6 hub genes were found to be VEGFA, KRAS, CDH2, ZEB2, TWIST1, and NTRK2 (Fig. 4B). We next constructed a circRNA-miRNA-hub gene sub-network. After excluding modules with inconsistent expression of circRNA and mRNA, 2 modules remained (hsa_circ_0000519/hsa-mir-204/CDH2 and hsa_circ_0000375/hsa-mir-3678/TWIST1) (Fig. 5).

Figure 4

Identification of hub genes from the PPI network. A. PPI network of 102 genes, consisting of 49 nodes and 55 edges. B. PPI network of 6 hub genes extracted from A. The node color changes gradually from yellow to red in ascending order according to the log2 (fold change) of genes.

Figure 5

CircRNA-miRNA-hub gene network. The network consists of 2 circRNAs, 2 miRNAs, and 3 hub genes. Vs stands for the circRNAs, while diamond is indicative of the miRNAs, rounded rectangles suggest the mRNAs, and grey edges represent the circRNA-miRNA-mRNA interactions. The red nodes indicate upregulated expression, whereas the green nodes indicate downregulated expression.

Functional enrichment analysis of DEmRNA

To identify the function of the 6 hub genes in the ceRNA network, GO enrichment KEGG pathway analysis was performed. A total of 507 GO terms, together with 15 KEGG pathways, were enriched with significant differences. The significant GO terms were “axonogenesis” in the biological process; “extrinsic component of plasma membrane” in the cellular components; and “neurotrophin binding” in the molecular function (P < .05). The top 5 GO terms for each category are indicated in Table 2. In addition, the top 10 KEGG analysis pathways included “Proteoglycans in cancer,” “Bladder cancer,” “Ras signaling pathway,” “VEGF signaling pathway,” “MAPK signaling pathway,” “MicroRNAs in cancer,” “Renal cell carcinoma,” “Pancreatic cancer,” and “EGFR tyrosine kinase inhibitor resistance” (Fig. 6).

Table 2

The top 5 GO terms enriched by DEmRNA involved in the ceRNA network.

Figure 6

Enrichment of the top 10 KEGG pathways of differentially expressed mRNAs in the ceRNA network. KEGG = Kyoto Encyclopedia of Gene and Genomes.

Identifying prognostic RNAs in patients

The miRNAs (hsa-mir-204 and hsa-mir-3678) and mRNAs (CDH2 and TWIST1) of the circRNA-miRNA-hub gene network were analyzed through the Kaplan–Meier method. High has-miR-204 expression was found in the Kaplan–Meier Plotter analysis and was attributed to a favorable OS (HR = 0.57, 95% CI: 0.46–0.7, P = 5.5e–08) (Fig. 7A). High has-miR-3678 was associated with poor OS (HR = 1.96, 95% CI: 1.42–2.73, P = 3.9e–05) (Fig. 7B). In addition, high expression of CDH2 was associated with poor OS (HR = 1.41, 95% CI: 1.14–1.75, P = .0016) and RFS (HR = 1.15, 95% CI: 1.03–1.28, P = .011) (Fig. 8A and B). However, no prognostic value of TWIST1 was found in OS and RFS (Fig. 8C and D).

Figure 7

The relationship between the expression of miRNAs and the overall survival of breast cancer patients. A, has-miR-204; B, has-miR-3678. The Kaplan–Meier curves were computed by utilizing the Kaplan–Meier Plotter web tool.

Figure 8

The relationship between the expression of mRNAs and the overall and relapse-free survival of breast cancer patients. CDH2: overall survival (A), relapse-free survival (B); TWIST1: overall survival (C), relapse-free survival (D). The Kaplan–Meier curves were computed by utilizing the Kaplan–Meier Plotter web tool.

Discussion

CircRNAs have been found to exhibit pivotal regulatory functions in carcinogenesis and cancer progression.[ Since circRNAs lack 5′ or 3′ polarities or polyadenylated tails, they are stable.[ In addition, circRNAs are abundant in eukaryotic cells and show a high degree of conservation, along with structural stability and a certain degree of organization, timing, and disease-specific activity.[ Based on these features, circRNAs have potential as biomarkers for cancers.[ Studies have confirmed dysregulated circRNAs in various types of cancer, and the circRNA-related ceRNA regulation network plays a crucial role in the pathogenesis and progression of tumors.[ However, the exact role of the circRNA-related ceRNA network in breast cancer development is still largely elusive. Therefore, comprehensive investigation of the impact of the ceRNA network impact on breast cancer is of critical importance. In this study, we first identified differentially expressed circRNAs, miRNAs, and mRNAs in breast cancer from public databases to provide a circRNA-miRNA-mRNA regulatory network.

Numerous studies have shown that circRNAs display dysregulated expression in breast cancer and are linked to pathogenesis and prognosis and are considered tumor-related biomarkers.[ Xu et al[ found that hsa_circ_001569 was upregulated in both breast cancer tissues and cell lines. Patients with increased hsa_circ_001569 possessed lymph-node metastasis, advanced clinical stage, and worse OS. Knockdown of hsa_circ_001569 reduced breast cancer cell growth and metastatic potential and activated PI3K-AKT signaling. Increasing evidence has shown that circRNA acts as a miRNA sponge during breast cancer development. For example, hsa_circ_001783 promoted the progression of breast cancer via sponging miR-200c-3p.[ In our study, we identified 2 key DEcircRNAs (hsa_circ_0000519 and hsa_circ_0000375) in the circRNA-miRNA-hub gene network. However, these 2 circRNAs have been previously reported in breast cancer.

MicroRNA (miRNA) is a highly conserved, regulatory single-stranded small RNA in the body. Although it does not directly encode proteins, it promotes the degradation of mRNA and inhibits protein translation, which in turn mediates the post-transcriptional regulation of genes.[ Abnormal expression of miRNAs has been shown to be involved in the development and progression of multiple tumors. MiRNAs, such as miR-181a, miR-373, and miR-490-3p, act as oncogenes that promote tumor cell proliferation, invasion, and metastasis.[ MiRNAs, such as miR-195-5p and miR-142-3p, can also act as tumor suppressors to promote tumor cell apoptosis and inhibit neovascularization.[ In the present study, we identified 2 miRNAs involved in the circRNA-miRNA-hub gene sub-network, including miR-204 and miR-3678. Of these 2 miRNAs, miR-204 has been reported in breast cancer.[ Hong et al[ indicated that miRNA-204-5p acts as a tumor suppressor to inhibit tumor metastasis and immune cell reprogramming in breast cancer by modulating PI3K/Akt signaling. Our results indicated that the expression of miR-204 was downregulated in breast cancer and associated with worse OS and RFS, consistent with previous studies.[ In addition, high miR-204 expression was associated with worse OS and RFS.

To further identify the key circRNAs participating in the regulatory network, we established the PPI network and selected 6 hub genes (VEGFA, KRAS, CDH2, ZEB2, TWIST1, and NTRK2). After excluding modules with inconsistent expression of circRNA and mRNA, the hsa_circ_0000519/hsa-mir-204/CDH2 and hsa_circ_0000375/hsa-mir-3678/TWIST1 axes were left. The CDH2 and TWIST1 hub genes play critical roles in the carcinogenesis and development of breast cancer. Our results showed that high expression of CDH2 was associated with poor OS and RFS. Lee et al[ found that miR-708-3p inhibits breast cancer cell epithelial-to-mesenchymal transition by targeting the CDH2 EMT activator. Strati et al[ indicated that elevated TWIST1 expression was associated with poor OS in EpCAM+ circulating tumor cells of early stage breast cancer patients. However, our study did not find prognostic value of TWIST1 in OS or RFS. Functional enrichment analyses indicated that these DEmRNAs have a significant effect on tumor-associated biological functions. Among the 15 pathways, “Ras signaling pathway,” “VEGF signaling pathway,” “MAPK signaling pathway,” “PI3K-Akt signaling pathway,” “EGFR tyrosine kinase inhibitor resistance,” “AGE-RAGE signaling pathway in diabetic complications,” and “Relaxin signaling pathway” are associated with the progression of breast cancer.[

Our study presents several limitations. First, the number of samples is not very large. Second, the conclusions of our study are only based on the current tools and databases. Third, the prognostic value of these DEcircRNAs in breast cancer has not been evaluated. In future studies, we will collect more clinical samples to validate our findings and further explore the function of these DEcircRNAs using in vitro and in vivo experiments.

Conclusions

We successfully established a ceRNA network to describe the possible mechanisms of breast cancer, which can shed light on the circRNA-related ceRNA network in breast cancer. The current study shows that hsa_circ_0000519 may play important roles in breast cancer. These findings may provide potential biomarkers or therapeutic targets for breast cancer patients.

Author contributions

Conceptualization: Chun-hua Zhao, Rui Qu.

Data curation: Chun-hua Zhao, Le Qu, Rui Qu.

Formal analysis: Chun-hua Zhao, Rui Qu.

Funding acquisition: Rui Qu.

Investigation: Chun-hua Zhao, Rui Qu.

Methodology: Chun-hua Zhao, Rui Qu.

Project administration: Chun-hua Zhao, Le Qu, Hui Zhang, Rui Qu.

Resources: Le Qu, Hui Zhang, Rui Qu.

Software: Le Qu, Hui Zhang, Rui Qu.

Supervision: Chun-hua Zhao, Rui Qu.

Validation: Chun-hua Zhao, Rui Qu.

Visualization: Chun-hua Zhao, Rui Qu.

Writing – original draft: Chun-hua Zhao, Le Qu, Hui Zhang, Rui Qu.

Writing – review & editing: Chun-hua Zhao, Le Qu, Hui Zhang, Rui Qu.