Construction of a circRNA-Related ceRNA Prognostic Regulatory Network in Breast Cancer

Introduction

Breast cancer is the second most common cancer and is the leading cause of cancer-related death among females worldwide with over 2 million newly diagnosed cases and more than 60 thousand deaths every year.¹ Despite advances in treatment, the mortality rate of breast cancer remains high, mainly due to its high frequency of metastasis, chemotherapy resistance, and recurrence. Some preventable causes of breast cancer have been identified, including aging, nulliparity, early menarche, late menopause, family history, and use of fertility drugs or oral contraceptives.^2–4 Recently, genetic and epigenetic factors have garnered increasing attention and have exhibited excellent application prospects in the diagnosis and treatment of breast cancer.

Noncoding RNA (ncRNA) is a functional RNA molecule that is transcribed from DNA but is not translated into proteins. Epigenetic related ncRNAs include tRNAs (transfer RNAs), rRNAs (ribosomal RNAs), miRNAs (microRNAs), siRNAs (small interfering RNAs), lncRNAs (long ncRNAs), circular RNAs (circRNAs), and piRNAs (Piwi-interacting RNAs). LncRNAs are characterized as being longer than 200 bp, while siRNAs, snRNAs, miRNAs, and piRNAs are shorter than 200 bp.⁵ Recently, ncRNAs have been shown to be involved in the pathogenesis and progression of various diseases through the mechanism of ceRNA function.^6–8 For instance, the lncRNA MYOSLID may function as a ceRNA to regulate MCL-1 expression by sponging miR-29c-3p in gastric cancer,⁹ and ciRS-7 could bind miR-7, serving as a therapeutic target by reducing EGFR-RAF1 activity in colorectal cancer.¹⁰

CircRNAs are a novel type of endogenous ncRNA that were considered an outcome of “splicing error” in the 1970s.^11,12 Advances in high-throughput sequencing technology have facilitated the understanding of circRNA function. CircRNAs can be classified as exotic, intronic, exon-intron, or intergenic.^11,13 Compared to corresponding linear RNAs, circRNAs are characterized by their covalently closed loop structure without a 5ʹ cap or a 3ʹ poly tail, which makes them resistant to RNase.¹⁴ Accumulating evidence has demonstrated specific expression patterns of circRNAs in cells, tissues, or different stages of diseases, presenting promising diagnostic and prognostic implications.^15,16 In terms of molecular mechanisms, circRNAs may function as competitive endogenous RNAs, they may bind proteins, and they may encode proteins,^17,18 participating in the occurrence and development of cardiovascular diseases, neurological disorders, and cancers.^19–21 However, the circRNA-related ceRNA network motifs and the underlying prognostic values of these RNAs in breast cancer have not been thoroughly elucidated.

In the current study, we combined our previously constructed circRNA database with the Gene Expression Omnibus (GEO) database and The Cancer Genome Atlas (TCGA) database to explore the circRNA-miRNA-mRNA axis in breast cancer.²² The GEO database and TCGA databases provide a large number of publicly available RNA-seq data regarding various cancers, including breast cancer.²³ Additionally, we attempted to construct a circRNA-related ceRNA regulatory network based on the dysregulated expression profiles in breast cancer, and we sought to further explore its pathogenesis and progression and to find the underlying prognostic and therapeutic biomarkers of breast cancer.

Materials and Methods

Patients and Data Collection

First, we recruited four patients with breast cancer who had undergone a mastectomy in the First People’s Hospital of Yixing City, China, in March 2016. We constructed circRNA expression profiles for breast cancer; the profiles were derived from four pairs of cancer tissues and adjacent normal-appearing tissues using Arraystar human circRNA microarray V2 (Catalog No: AS-CR-H-V2.0).²⁴ The circRNA expression profiles were recruited from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) with the filter criteria of 1) the expression dataset was acquired from tumor and normal tissues in patients with breast cancer; 2) the detection method was microarray analysis; and 3) the sample size was at least three. Finally, the dataset GSE101123, consisting of eight breast cancer and three normal tissues, was selected. The level-three RNA-seq raw count data of miRNA (containing 1090 breast cancer and 103 normal tissues) and mRNA (including 1096 breast cancer and 112 normal tissues) were extracted from the TCGA data portal (https://portal.gdc.cancer.gov/, access date: August 14, 2019), except for the expression data of males. Additionally, we downloaded the corresponding clinical characteristics from TCGA, including age, survival time, survival status, clinical stage, and T, N, and M stages. The flow chart of the analysis is illustrated in Figure 1.

Figure 1 Flowchart of circRNA-related ceRNA regulatory network analysis in breast cancer.

Identification of Differentially Expressed circRNAs, miRNAs and mRNAs

First, we extracted and sorted the expression matrixes by Perl (https://www.perl.org/get.html) and identified the differential expression using R3.6.0 (https://www.r-project.org/). For differentially expressed circRNA (DEcircRNA) analysis, we first integrated our previous circRNA dataset with GSE101123 by batch normalization using the “sva” package and then applied the “limma” package to screen DEcircRNAs with a filter criterion of |Fold Change| > 1 and FDR < 0.05. The filter criteria for differentially expressed miRNAs (DEmiRNAs) and mRNAs (DEmRNAs) were |log2(fold-change)| > 2 and FDR < 0.05 using the “edgeR” package. We applied the “pheatmap” package to construct the heatmap of DEcircRNAs, DEmiRNAs, and DEmRNAs. Afterward, we established the coexpression network for DEmRNAs in breast cancer using a weighted gene correlation network analysis (WGCNA) package. First, we used goodSampleGenes and meanFPKM =0.5 to remove the unqualified genes. Then, we picked the best soft-thresholding power by the pickSoftThreshold function. Adjacency was converted into a TOM, and genes with high correlation were grouped into different modules according to TOM-based dissimilarity with minModuleSize =30. Ultimately, the DEmRNAs involved in the highest correlated modules were candidates for further investigation, and we performed logarithmic transformation for homogenization.

Construction of the circRNA-miRNA-mRNA Regulatory Network

We employed the cancer-specific circRNA database (CSCD, http://gb.whu.edu.cn/CSCD/) and the circular RNA interactome (CircInteractome, https://circinteractome.nia.nih.gov/) database to predict circRNA-miRNA interactions. The target miRNAs were further identified by overlapping with DEmiRNAs obtained from the TCGA database. Pairs of miRNA-mRNA interactions were predicted by three bioinformatic databases, namely, miRDB (https://mirdb.org/download.html), miRTarBase (http://mirtarbase.mbc.nctu.edu.tw/php/index.php) and TargetScan (http://www.targetscan.org/vert_72/). Only the interactions that were consistently predicted by the three databases were selected for further analysis. Then, we screened the target mRNAs by intersecting with the DEmRNAs involved in the key modules based on WGCNA. Finally, based on miRNA-mRNA pairs and circRNA-miRNA interactions, we established the circRNA-miRNA-mRNA regulatory network and visualized it with Cytoscape 3.7.1 software (https://cytoscape.org/).

Identification of Prognostic Signatures in Breast Cancer

To investigate the prognostic value of targeted mRNAs involved in the circRNA-related ceRNA network of breast cancer, we constructed a prognostic model. Briefly, we identified the prognostic signatures based on Cox regression analysis using the “survival” package and selected the key signatures according to the best Akaike information criterion (AIC) using the “MASS” package. These signatures were subsequently applied to construct a prognostic model and generate a risk score based on the expression of each signature and its corresponding coefficients. We also performed a receiver operating characteristic curve (ROC) to evaluate the risk score. Patients with risk scores lower than the median values were classified into a low-risk group, and the rest were classified into a high-risk group.

Functional Analysis of Prognostic Signatures in the circRNA-Related Network

To predict the potential functions of prognostic signatures involved in the circRNA-related ceRNA network, we performed Gene Set Enrichment Analysis (GSEA, https://www.gsea-msigdb.org/gsea/index.jsp).²⁵ For the phenotypic file, we defined patients with targeted mRNA expression greater than the median level as the “high-risk” subgroup and defined the rest as the “low-risk” subgroup. We applied the C2 (c2.cp.kegg.v7.1.symbols.gmt) subcollection as reference gene sets. Thresholds for significance were determined by 1000 permutation analyses, and FDR < 0.05 was considered to be significant.

Ethics Statement

This study was approved by the Ethics Committee of Nanjing Medical University. After informed consent was obtained from all participants, questionnaires were used to collect demographic data. This study was conducted in accordance with the Declaration of Helsinki. Data downloaded from the public database followed the data access policy.

Results

Differentially Expressed circRNAs, miRNAs, and mRNAs in Breast Cancer

Compared with normal breast tissues, a total of 59 DEcircRNAs were identified from circRNA microarray analyses of breast-cancer tissues, including 54 upregulated circRNAs and 5 downregulated circRNAs (Figure 2A). The basic characteristics of the top 4 DEcircRNAs (hsa_circ_0062682, hsa_circ_0007798, hsa_circ_0092342, and hsa_circ_0065173) are described in Table 1. A total of 98 DEmiRNAs were identified, including 78 upregulated miRNAs and 20 downregulated miRNAs (Figure 2B). Among the DEmiRNAs, miR-133b (log2(fold-change) = −6.96, FDR < 0.05) was the most downregulated, and miR-122-5p (log2(fold-change) = 7.79, FDR < 0.05) was the most upregulated. Meanwhile, 3966 DEmRNAs, consisting of 2906 upregulated mRNAs and 1060 downregulated mRNAs, were detected between breast-cancer tissues and normal tissues (Figure 2C), of which UCN3 (log2(fold-change) = 9.09, FDR < 0.05) and MYH2 (log2(fold-change) = −8.59, FDR < 0.05) were ranked as the most upregulated and downregulated mRNAs, respectively.

Table 1 Characteristics of the Top Four DEcircRNAs in Breast Cancer

Figure 2 Differentially expressed genes in breast cancer. (A) Heatmap and volcano of DEcircRNAs identified from previous circRNA datasets and GEO databases using the limma package (|log2(fold-change) >1| and FDR < 0.05); (B) Heatmap and volcano of DEmiRNAs identified from TCGA database using the edgeR package (|log2(fold-change) >2| and FDR < 0.05); (C) Heatmap and volcano of DEmRNAs identified from TCGA database using the edgeR package (|log2(fold-change) > 2| and FDR < 0.05).

Identification of the Coexpression Network for DEmRNAs in Breast Cancer

To identify the key modules of mRNAs in breast cancer, we performed WGCNA on the DEmRNAs. In brief, β =4 was selected as the best soft threshold with a scale-free R² = 0.96. Subsequently, the highly correlated modules were merged with 75% relevance, and a total of 11 coexpression modules were generated, with the number of mRNAs contained in these modules ranging from 44 to 1400 (Figure 3A). As shown in Figure 3B, the most significant module associated with breast cancer and the normal trait was the turquoise module, which was composed of 1400 mRNAs that were candidates for further analysis (P = 2×10⁻²⁶²). The scatterplots of Gene Significance vs Module Membership in the turquoise module suggested the highest correlation (Figure 3C).

Figure 3 Weighted gene correlation network analysis (WGCNA) of mRNAs in breast cancer. (A) Identification of a coexpression module of DEmRNAs in breast cancer. (B) The correlation between the DEmRNA module and clinical traits. (C) Scatter plot of module eigengenes in the turquoise modules. P < 0.05. Abbreviation: DEmRNAs, differentially expressed mRNAs.

Construction of the ceRNA Regulatory Network in Breast Cancer

To clarify the circRNA-related ceRNA regulatory mechanism, a circRNA-miRNA-mRNA network was constructed. We used the CSCD database for the miRNA response element (MRE) prediction of DEcircRNAs, and the CircInteractome database was subsequently applied for those not included in the CSCD. After intersecting with the DEmiRNAs identified from the TCGA database, a total of 46 miRNAs were identified. Next, we acquired 1761 predicted mRNAs of these 46 miRNAs by overlapping the three databases (miRDB, miRTarBase, and TargetScan). A total of 30 predicted mRNAs remained after overlapping with DEmRNAs involved in key modules of DEmRNAs by WGCNA. Finally, a total of 38 circRNA-miRNA pairs, including 29 circRNAs and 18 miRNAs, and 38 miRNA-mRNA pairs, including 18 miRNAs and 30 mRNAs, remained for further ceRNA network construction. The ceRNA network based on 76 intersections is illustrated in Figure 4A.

Figure 4 CircRNA-miRNA-mRNA ceRNA regulatory network in breast cancer. (A) The ceRNA regulatory network in breast cancer. (B) The prognostic ceRNA regulatory network in breast cancer.

Exploration of the Prognostic Model and Corresponding ceRNA Network

To estimate the prognostic values of mRNAs involved in the circRNA-related ceRNA regulatory network, we attempted to construct a prognostic model. A total of 30 mRNAs involved in the circRNA-related ceRNA network were applied to establish the prognostic model. Then, 7 signatures were identified as optimally prognostic signatures based on stepAIC function, including MMD, SLC29A4, CREB5, FOS, ANKRD29, MYOCD and PIGR. Kaplan–Meier analysis of these 7 prognostic signatures for patients with breast cancer is shown in Supplemental Figure 1. Lower expression of ANKRD29 and PIGR was associated with poor prognosis in patients with breast cancer (P = 0.0371 and P = 0.00632, respectively).

After that step, a risk score for the prediction of prognosis was generated as follows: Risk score=Expression_MMD*0.21041-Expression_SLC29A4*0.11003-Expression_CREB5*0.14338 -Expression_FOS*0.07729-Expression_ANKRD29*0.12044+Expression_MYOCD*0.16679-Expression_PIGR*0.09200. We found that the risk score can accurately predict the prognosis of patients with an area under the curve (AUC) = 0.833 (Figure 5A). Next, we classified patients into high- and low-risk groups based on the median risk score. As shown in Figure 5B, the overall survival rate of patients in the high-risk group was lower than that of patients in the low-risk group (P = 1.356×10⁻⁶). Based on the univariate Cox regression analysis, the age (HR:1.036, 95% CI: 1.021–1.051), clinical-stage (HR:2.768, 95% CI: 1.932–3.964; stage III-IV vs stage I-II), T stage (HR:1.925, 95% CI: 1.291–2.869; T3-4 vs T1-2), M stage (HR:5.907, 95% CI: 3.244–10.757; M1 vs M0), N stage (HR:2.162, 95% CI: 1.486–3.147; N1-3 vs N0) and risk score (HR:2.469, 95% CI: 1.680–3.628) were related to the overall survival (Figure 5C). As shown in Figure 5D, the multivariate Cox regression analysis indicated that the risk score was an independent risk factor for overall survival (HR: 2.305, 95% CI: 1.533–3.464). In addition, we observed that 7 signatures were involved in 20 circRNA-related ceRNA networks (Figure 4B).

Figure 5 Construction of 7 signature-based prognostic models. (A) The ROC curve to evaluate the risk score. (B) Kaplan–Meier analysis for patients with high and low risk scores for breast cancer. (C) Univariate Cox regression analysis. (D) Multivariable Cox regression analysis. P < 0.05. Abbreviation: ROC, receiver operating characteristic curve.

GSEA Analysis for Prognostic Signatures in the ceRNA Regulatory Network

We performed GSEA analysis on 1096 patients from high-risk and low-risk groups to identify the cancer-related pathways involved in the 7 prognostic signatures. A total of 27 pathways were enriched in patients with high expression of MMD, 1 pathway for SLC29A4, 4 pathways for CREB5, 21 pathways for FOS, 14 pathways for ANKRD29, 30 pathways for MYOCD and 9 pathways for PIGR. As shown in Figure 6, multiple cancer-related pathways were enriched, including the TGF-β pathway (NES = 2.35, FDR < 0.001) for MMD, ECM receptor interaction for CREB5 (NES = 2.37, FDR < 0.001), ECM receptor interaction for FOS (NES = 2.09, FDR<0.001), cell cycle for ANKRD29 (NES = −2.27, FDR<0.001), cell cycle for MYOCD (NES = −2.18, FDR < 0.001) and JAK STAT signaling pathway (NES = −1.90, FDR = 0.04) for PIGR.

Figure 6 GSEA analysis of these 7 signatures in breast cancer. (A) GSEA analysis of MMD; (B) GSEA analysis of SLC29A4; (C) GSEA analysis of CREB5; (D) GSEA analysis of FOS; (E) GSEA analysis of ANKRD29; (F) GSEA analysis of MYOCD; (G) GSEA analysis of PIGR. (FDR < 0.05).

Discussion

Breast cancer is one of the most common malignant tumors among females in the world. It is imperative to clarify the molecular mechanisms of carcinogenesis, explore prognostic factors, and search for therapeutic biomarkers of breast cancer. NcRNAs are a novel type of RNA that is primarily expressed in eukaryotes and plays an essential role in maintaining regular physiological activity.^26–30 Although the ceRNA regulatory network has been widely verified, the circRNA-related ceRNA network has rarely been studied. In the current study, we established a circRNA-related ceRNA regulatory network and prognostic model in breast cancer by bioinformatic analysis. Expression datasets of DEcircRNAs, DEmiRNAs, and DEmRNAs between breast-cancer tissues and adjacent normal mammary glands were identified based on our previous circRNA dataset, the GEO database and TCGA database.^31,32 Then, we constructed a circRNA-miRNA-mRNA regulatory network to investigate the underlying molecular mechanisms of DEcircRNAs. We constructed a prognostic model derived from 7 signatures in the ceRNA network, and functional analysis suggested that the 7 signatures were significantly enriched in cancer-related pathways, thereby demonstrating promise as prognostic and therapeutic biomarkers of breast cancer.

With the rapid development of high-throughput sequencing technology, the biogenesis and function of circRNAs have been widely explored. Accumulating evidence has confirmed the capacity of circRNAs to serve as a biomarker in the diagnosis or prognosis of cancer due to its covalently closed loop structure. For instance, a decrease in circMTO1 can predict poor survival among patients with hepatocellular carcinoma.³³ Wang et al reported that hsa_circ_0077837 and hsa_circ_0001821 could serve as potential biomarkers for lung cancer, whereas hsa_circ_0001073 and hsa_circ_0001495 showed value as diagnostic molecules.³⁴ Circ-KIAA1244, a gastric cancer-derived circRNA, was reported to be a novel circulating molecular biomarker of gastric cancer.³⁵ Xu et al reported that hsa_circ_0005230 could mechanically sponge miR-1238 and miR-1299 to serve its oncogenic functions in cholangiocarcinoma.³⁶ In the current study, we identified 59 DEcircRNAs in breast cancer, including hsa_circ_0062682, hsa_circ_0007798, hsa_circ_0092342, and hsa_circ_0065173. Among these DEcircRNAs, hsa_circ_0065173 has been previously found to be altered in breast cancer.³⁷

Several molecular functions of circRNAs have been identified, including sponging microRNA, binding to RNA-binding proteins and serving protein translation functions.^38,39 CircRNAs are enriched with miRNA binding sites, which causes them to act as miRNA sponges and then participate in the pathogenesis and development of various diseases. Zhang et al reported that the downregulation of circRNA_100269 in gastric cancer may suppress cell line proliferation (in AGS and MKN28 cells) by targeting miR-630.⁴⁰ CircRNA ZFR could stimulate the progression of non-small-cell lung cancer by acting as a miR-101-3p sponge to enhance CUL4B expression.⁴¹ In this study, we successfully constructed a circRNA-miRNA-mRNA axis, including 29 DEcircRNAs, 18 DEmiRNAs, and 30 DEmRNAs, which might help to elucidate the circRNA-related regulatory mechanisms in breast cancer. In addition, we successfully constructed a molecular prognostic model originating from 7 signatures involved in the circRNA-related ceRNA regulatory network (MMD, SLC29A4, CREB5, FOS, ANKRD29, MYOCD and PIGR), providing new insight into postoperative management for breast cancer. Meanwhile, functional analysis indicated that these 7 signatures were significantly enriched in cancer-related pathways, including the TGF-β pathway, the focal adhesion pathway, ECM-receptor interaction, the cell cycle pathway and the JAK-STAT signaling pathway. We confirmed 7 signatures in 20 ceRNA networks, which might help elucidate the mechanisms of circRNA in breast cancer.

In this study, we combined the publicly available TCGA and GEO databases with our previous circRNA dataset for bioinformatic analysis, identified 59 DEcircRNAs, 98 DEmiRNAs, and 3966 DEmRNAs and constructed a breast cancer-related ceRNA regulatory network, including 38 circRNA-miRNA pairs and 38 miRNA-mRNA pair intersections. We also successfully generated a prognostic model based on 7 signatures and 20 prognostic ceRNA axes. Multiple cancer-related pathways were involved in the ceRNA network, providing new insight into the underlying pathogenesis of breast cancer and suggesting therapeutic targets for patients with this disease. The results of this study may help to deepen the understanding of circRNA-related regulatory mechanisms. Moreover, we generated a prognostic model that provided new insight into postoperative management for breast cancer.

However, it should be noted that our findings are based on a microarray analysis of a small sample followed by bioinformatic analysis. Further in-depth studies should be conducted to verify the underlying ceRNA mechanism in breast cancer.