This article is mentioned in:

Introduction

Breast cancer (BC) is the most common type of cancer among females and accounts for ~16% of all cancers (1). BC is a heterogeneous malignant tumor derived from breast tissue and is the main cause of cancer-associated mortality in females worldwide (2,3). Surgery and adjuvant therapy can effectively improve the prognosis of patients with BC; however, metastasis and recurrence can still occur (4,5). The development of BC is a complex process, involving genetic and epigenetic changes (6). Due to the heterogeneity of BC, clinical outcomes vary. The predictive value of current diagnostic and prognostic factors for BC is limited. Thus, it is essential to identify effective diagnostic markers, therapeutic targets and prognostic factors. Recently, competing endogenous RNAs (ceRNAs) have been investigated in the process of tumorigenesis and cancer development (7,8). Salmena et al (7) proposed the ceRNA hypothesis in 2011, in which a complex post-transcriptional control network involving all types of RNA transcripts suppresses the functions of microRNAs (miRNAs/miRs) via miRNA response elements.

Long non-coding RNAs (lncRNAs) are a class of noncoding RNAs with a length >200 nucleotides (9); a lncRNA may function as an oncogene or a tumor suppressor to regulate the expression of target genes at various levels (10,11). Increasing evidence has demonstrated that lncRNAs serve an important role in several cellular processes (12,13). The predictive value of lncRNAs has been reported in multiple types of cancer, including lung cancer (14), gallbladder cancer (15), ovarian cancer (16), glioblastoma (17), gastric cancer (18), pancreatic cancer (19) and colorectal cancer (20). Abnormal expression of certain lncRNAs has been observed in BC, and it has been associated with diagnosis and prognosis of the disease (21). However, due to the limited number of tissue samples employed, the results of these studies in BC are inconsistent (22,23).

In the present study, the expression profiles of mRNAs, miRNAs and lncRNAs in BC and non-cancer tissues were obtained from The Cancer Genome Atlas (TCGA). An lncRNA-ceRNA network associated with BC was successfully constructed. In addition, functional enrichment analysis was conducted to investigate the potential mechanism. Survival analysis was also performed for the differentially expressed genes (DEGs) and their association with BC prognosis was explored. Differentially expressed RNAs (DERNAs) in the ceRNA network were analyzed to further investigate the development and prognosis of BC.

Materials and methods

Data acquisition and processing

RNA sequencing (RNA-Seq) and the corresponding clinical data of BC patients were obtained from the TCGA data portal (https://tcga-data.nci.nih.gov/tcga/). The lncRNA, mRNA and miRNA sequence data were derived from Illumina HiSeq RNA Seq and Illumina HiSeq miRNA Seq platforms. A total of 1,109 BC tissues and 113 normal breast tissues were included in the present study. The inclusion criteria were as follows: i) Gender was female; and ii) histopathology was confirmed as ductal or lobular disease. The exclusion criteria were as follows: i) Gender was male; and ii 2) histopathology was not confirmed as ductal and lobular disease. The present study was conducted in accordance with the publication guidelines provided by the TCGA (http://cancergenome.nih.gov/publications/publicationguidelines). Therefore, further approval from an ethics committee was not required.

Identification of DERNAs

Differentially expressed lncRNA (DElncRNA) and differentially expressed mRNA (DEmRNA) were defined and encoded based on the annotations from the Ensembl database (http://www.ensembl.org/index.html). To further analyze the data, the edgeR package in the R language 3.53 was used (http://bioconductor.org/packages/edgeR/), and the DElncRNA, DEmRNA and differentially expressed miRNA (DEmiRNA) were identified. A |log2 fold change (FC)| >1.5 and false discovery rate (FDR) adjusted to a P-value <0.01 were set as the thresholds. Volcano plots and heatmaps of the DERNAs were produced by using the gplots and heatmap packages in the R platform, respectively.

Construction of the ceRNA network

The miRcode (http://www.mircode.org/) database (24) was used to predict lncRNA-miRNA interactions, and the StarBase v2.0 (25) database (http://starbase.sysu.edu.cn/) was used to decode the miRNA sequences. In addition, miRNA-targeted mRNAs were retrieved from the miRDB (26), miRTarBase (27), and TargetScan (28) databases. Only the mRNAs recognized by all three databases were considered as candidate mRNAs and were then intersected with the DEmRNAs to determine the DEmRNAs that were targeted by the DEmiRNAs. Based on the DEmiRNA-DElncRNA and DEmiRNA-DEmRNA interactions, a coexpression network of DEGs was constructed, which was visualized using Cytoscape 3.6.1 (29). In addition, the cytoHubba plugin was used to evaluate the top 15 genes in the network based on closeness.

Functional enrichment analysis

To understand the underlying biological mechanisms of the DEmRNAs in the ceRNA crosstalk network, Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were conducted using the DAVID (database for annotation, visualization, and integrated discovery; http://david.ncifcrf.gov/) online tool (30) and cluster profiler, which is an R package for functional classification and enrichment of gene clusters using hypergeometric distribution. The GO plot package of R software was used to display the results of the GO and KEGG analyses. GO and KEGG enrichment analyses were based on the threshold of P<0.05.

Survival analysis

The high expression and low expression groups were classified according to the median. Survival analysis was performed for the DERNAs in the ceRNA network using the ‘survival’ package in R to assess their prognostic value in patients with BC. Survival curves were generated using the Kaplan-Meier method, and the log-rank test was used to compare the differences between the groups. P<0.05 was considered as the threshold for significance in all analyses.

Results

DEmRNAs, DElncRNAs and DEmiRNAs in BC

Significant DEGs were identified in 1,109 BC and 113 adjacent non-cancer breast tissues using the ‘edgeR’ package in R software. A total of 3,198 DEmRNAs (1,996 upregulated and 1,202 downregulated), 1,043 DElncRNAs (814 upregulated and 229 downregulated), and 150 DEmiRNAs (114 upregulated and 36 downregulated) were identified with thresholds of |log2 fold change (FC)| >1.5 and adjusted P<0.01. The distribution of DEGs based on the two parameters, -log false discovery rate and logFC, was presented as a volcano plot (Fig. 1A). The heatmap demonstrated the expression data of all DEGs (Fig. 1B).

Figure 1. Gene expression analysis of BC data from The Cancer Genome Atlas database. (A) Volcano plots of differentially expressed lncRNAs, miRNAs and mRNAs. The red points in the plots represent upregulated RNAs and the green points represent downregulated RNAs with statistical significance. (B) Heatmaps of differentially expressed RNAs between BC and non-tumor tissues. BC, breast cancer; lncRNA, long non-coding RNA; miRNA, microRNA.

Construction of the ceRNA network

To investigate how lncRNAs mediate mRNA expression by integrating with miRNAs in BC, a lncRNA-miRNA-mRNA (ceRNA) network was generated and visualized using Cytoscape v3.6.1. A total of 405 pairs of interacting lncRNAs and miRNAs were identified using the Perl program 5.30.0 (http://www.perl.org/) (31), based on the 1,043 DElncRNAs. Based on 24 miRNAs, target mRNAs were identified using the miRTarBase, miRDB and TargetScan databases. Target genes that were included in all three datasets were selected. Finally, there were 74 DEmRNAs in the ceRNA network (Fig. 2). In total, 97 lncRNA nodes, 24 miRNA nodes and 74 mRNA nodes were identified as having differential expression profiles in the ceRNA network (Fig. 3A). Based on the plug-in cytoHubba in Cytoscape, the hub network was obtained through the existing network analysis and closeness scores (Fig. 3B). The results included six miRNAs, namely hsa-miR-204, hsa-miR-145, hsa-miR-122, hsa-miR-155, hsa-miR-182 and hsa-miR-183. The lncRNAs involved included LINC00461, AGAP11, LINC00466, LINC00261, MAGI2-AS3, TCL6, ADAMTS9-AS2, AL589642.1 and DLX6-AS1.

Figure 2. Venn diagram of differentially expressed mRNAs involved in the competing endogenous RNA network. The red circle represents the differentially expressed genes, and the blue circle represents the target genes of miRNAs in the competing endogenous RNA network.

Figure 3. Schematic of the ceRNA network of lncRNA-miRNA-mRNA in breast cancer. (A) The complete network is shown, with diamonds indicating IncRNAs, rounded rectangles indicating miRNAs, and ellipses indicating mRNAs. The red and blue nodes represent upregulation and downregulation, respectively. (B) The hub closeness network involved in the ceRNA network. ceRNA, competing endogenous RNA; lncRNA, long non-coding RNA; miRNA, microRNA.

Functional analysis of DEmRNAs in the ceRNA network

Using GO and KEGG analysis, the biological functions of the 74 DEmRNAs were investigated. The results indicated that the DEmRNAs were enriched in two GO biological process categories (P<0.05; Fig. 4A). The main biological process GO term was ‘positive regulation of GTPase activity’. In addition, enriched GO terms based on the analysis of survival curves associated with KM-significant DEmRNAs were identified (Fig. 4B). The network included nine DEmRNAs, which were mainly enriched in ‘positive regulation of epithelial cell proliferation’.

Figure 4. GO enrichment analysis. (A) Enriched GO biological process terms of DEmRNAs involved in the ceRNA network. (B) Enriched GO biological process terms of DEmRNAs involved in the ceRNA network and significantly affecting prognosis by survival analysis. GO, gene ontology; DE, differentially expressed; ceRNA, competing endogenous RNA.

From KEGG analysis, four main enriched pathways were determined (Table I and Fig. 5), including the ‘mitogen-activated protein kinase (MAPK) signaling pathway’, the ‘Ras signaling pathway’, ‘prostate cancer’ and the ‘FoxO signaling pathway’. Among the four pathways, the ‘MAPK pathway’ and the ‘Ras pathway’ are associated with BC progression (3,30). In addition, ‘prostate cancer’ is linked to tumor-related pathways (32), whereas the ‘FoxO signaling pathway’ is associated with apoptosis, cell-cycle control and glucose metabolism (33).

Figure 5. Top 4 enriched Kyoto Encyclopedia of Genes and Genomes pathways of differentially expressed mRNAs involved in the competing endogenous RNA network. MAPK, mitogen-activated protein kinase; FoxO, forkhead box O.

Table I. Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis of the differentially enriched mRNAs.

Table I.

Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis of the differentially enriched mRNAs.

Pathway ID	Description	Genes involved	Count	P-value	Adjusted P-value (FDR)
hsa04010	MAPK signaling pathway	TGFBR2, KIT, FGF2, PDGFD, FOS, NTRK2, FGFR3, RASGRP1	8	<0.001	0.014
hsa05215	Prostate cancer	FOXO1, ERG, CREB5, PDGFD, CCNE2	5	<0.001	0.014
hsa04014	Ras signaling pathway	KIT, FGF2, RGL1, PDGFD, NTRK2, FGFR3, RASGRP1	7	<0.001	0.014
hsa04068	FoxO signaling pathway	TGFBR2, FOXO1, GABARAPL1, S1PR1, CCNB1	5	<0.001	0.032

Survival-associated DERNAs in the ceRNA network

KM analysis was performed to investigate the association between DERNAs (DElncRNAs, DEmRNAs and DEmiRNAs) in the ceRNA network and the prognosis of patients with BC, thereby obtaining a prognostic signature for each DERNA.

In total, 6 out of the total 97 DElncRNAs were significantly associated with the rate of overall survival. Among them, three DElncRNAs (INC AC112721.1, LINC00536 and MIR7-3HG) were negatively associated with OS, while the other three (ADAMTS9-AS1, AL356479.1 and LINC00466) were positively associated with OS (Fig. 6A).

Figure 6. Kaplan-Meier survival curves of differentially expressed (A) long non-coding RNAs, (B) mRNAs and (C) microRNAs significantly associated with overall survival in breast cancer. Table I. Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis of the differentially enriched mRNAs.

In addition, 9 out of the total 74 DEmRNAs were significantly associated with the rate of overall survival. Among them, four DEmRNAs [importin subunit α 2 (KPNA2), Rac GTPase-activating protein 1 (RACGAP1), Shc binding and spindle associated 1 (SHCBP1) and zinc finger protein 367 (ZNF367)] were negatively associated with OS, while the remaining five DEmRNAs [tropomyosin receptor kinase B (NTRK2), ORS1, prostaglandin-endoperoxidase synthase 2 (PTGS2), RAS guanyl-releasing protein 1 (RASGRP1) and secreted frizzled-related protein 1 (SFRP1)] were positively associated with OS (Fig. 6B).

Finally, 2 out of the total 24 DEmiRNAs were significantly associated with the rate of OS. The results demonstrated that hsa-miR-301b was negatively associated with OS, while hsa-miR-204 was positively associated with OS (Fig. 6C).

Discussion

BC is a common gynecological malignancy; the incidence of this disease is increasing (3). Despite the combined use of surgery, chemotherapy, radiotherapy and endocrine therapy, BC remains a major cause of cancer-associated mortality in females (34). To improve the clinical status of patients with BC, and develop additional methods for the diagnosis, treatment and prognosis of this disease, it is crucial to understand the regulatory mechanism underlying its occurrence and development. The study of prognostic markers associated with BC may aid early detection and individualized treatment.

Recently, research into the lncRNA-miRNA-mRNA network has attracted increasing attention. Non-coding (nc) RNAs include various forms of RNA, which are categorized into lncRNAs and short ncRNAs. As short-chain ncRNAs, miRNAs are ~22 nucleotides in length and bind to sequences with partial complementarity to RNA transcripts to act as negative regulators of target gene expression (35). lncRNAs are defined as transcripts >200 nucleotides in length with no protein-coding functions (36). Denzler et al (37) proposed the ‘competitive endogenous RNA’ hypothesis, suggesting that ceRNAs could regulate miRNAs. Salmena et al (7) reported that all types of RNAs can compete with each other for miRNAs, resulting in large-scale transcription crosstalk throughout the transcriptome.

Investigation into ceRNAs and ceRNA networks may provide insight into the processes underlying certain diseases and aid the development of novel therapeutic strategies. The imbalance in lncRNA expression is associated with tumorigenesis, metastasis and prognosis of various types of cancer (14,16,20,38), suggesting that lncRNAs may be potential biomarkers of cancer. In order to evaluate the systematic biological function of the DElncRNAs in the ceRNA network, the present study investigated the association between DElncRNAs and OS. To comprehensively understand how the lncRNA-associated ceRNA network affects BC, large-scale sequencing data of patients with BC were obtained for analysis from TCGA. DERNAs were detected and a lncRNA-miRNA-mRNA network was constructed based on the prediction of associated tumor-specific biological processes. In total, 97 lncRNA, 22 miRNA and 74 mRNA nodes in the ceRNA networks were identified to be differentially expressed. In addition, the association between DElncRNAs and OS, and the biological role of DElncRNAs in ceRNA networks were investigated. In the present study, nine DEmRNAs (KPNA2, NTRK2, ORS1, PTGS2, RACGAP1, RASGRP1, SFRP1, SHCBP1 and ZNF367), six DElncRNAs (INC AC112721.1, ADAMTS9-AS1, AL356479.1, LINC00466, LINC00536 and MIR7-3HG), and two DEmiRNAs (hsa-miR-204, hsa-miR-301b) were associated with OS; thus, these DERNAs may serve as potential prognostic biomarkers of OS in patients with BC.

The regulation of miRNAs is a critical process associated with tumor-related pathways involved in oncogenesis and tumor progression (39). It has been reported that a higher grade of clear cell renal cell carcinoma is correlated with a concomitant decrease in miR-204 expression (40,41). Ying et al (42) revealed that miR-204 suppresses the stem cell-associated phenotype, and the self-renewal and migration potentials of glioma cells by targeting SRY-box transcription factor 4 and the EPH receptor B2. Thus, loss of miR-204 expression enhances tumor migration. Shen et al (43) suggested that miR-204 regulates the biological behavior of BC cells by directly targeting forkhead box A1, acting as a tumor suppressor in BC. These reports are consistent with the current results, in which the increased expression of miR-204 was associated with improved overall survival in patients with BC. The current results also demonstrated that upregulated hsa-miR-301b was associated with poor survival in patients with BC. Fort et al (44) revealed that the miR-130b/301b cluster serves a protumorigenic role in prostate cancer.

In the present study, increased expression of ADAMTS9-AS1, AL356479.1 and LINC00466 was positively associated with prognosis, while the reduced expression of the remaining three lncRNAs (INC AC112721.1, LINC00536 and MIR7-3HG) was associated with improved prognosis. As an antisense lncRNA, ADAMTS9-AS1 serves crucial roles in the progression and prognosis of carcinoma (45). ADAMTS9-AS1 was determined to be associated with the prognosis of patients with bladder cancer (46,47). Xing et al (48) suggested that ADAMTS9-AS1 may be a marker for predicting the prognosis of colorectal carcinoma. Li et al (49) reported that ADAMTS9-AS1 may be indispensable in the development of the ectoderm and epithelial cells, and can be effectively employed for determining the risk of esophageal cancer development. In addition, it was indicated that ADAMTS9-AS1 serves as a risk-predicting lncRNA in ovarian cancer (50). Analysis of BC data revealed ADAMTS9-AS1 and LINC00536 to be prognostic markers for BC (50). These previous findings were consistent with the results of the present study. To the best of our knowledge, no studies have reported the roles of AC112721.1, AL356479.1, LINC00466 and MIR7-3HG in cancer.

Accumulating evidence has demonstrated that miRNAs serve an important role in cancer. In the ceRNA network, the present study reported that the upregulated expression of LINC00466 was associated with downregulation of hsa-miR-204 and hsa-miR-144; and MIR7-3HG was determined to downregulate hsa-miR-204 and hsa-miR-145. Several studies have reported that the levels of hsa-miR-204 were decreased in certain types of cancer, including gastric cancer (51) and lung cancer (52). Based on TCGA data of 1,077 BC tissue samples, Alhawaj et al (53) determined that miR-204-5p was markedly downregulated when compared with 104 normal tissue samples. Liu et al (54) demonstrated that trichostatin A and tamoxifen inhibited BC cell growth via miR-204 and estrogen receptor (ER) α-mediated suppression of the AKT/mTOR pathway, suggesting that miR-204 may be associated with tumor growth and with sensitivity to chemotherapeutic agents. The present study proposed the involvement of the LINC00466-hsa-miR-204-NTRK2 axis in BC. Thus, LINC00466 and MIR7-3HG may be crucial for the progression of tumors by interacting with hsa-miR-204. In addition, hsa-miR-144 and hsa-miR-145 exhibit inverse expression profiles in BC tissues (55–57).

Combining the results of the hub networks based on analysis with the plug-in cytoHubba in Cytoscape, the present study determined that hsa-miR-155 and AGAP11 had a high degree of closeness, and exhibited significance from GO analysis. Hsa-miR-155 is a multifunctional small RNA located in chromosome 21, and is involved in a variety of biological and pathological processes (58). Volinia et al (59) reported that miRNAs served complex roles in tumors; network analysis indicated that the majority of miRNAs were linked to hsa-miR-155 within transgenic mice with miR-155-induced acute lymphocytic leukemia (59). Chuang et al (60) identified that the increased expression of hsa-miR-155-5p was a poor prognostic marker in acute myeloid leukemia. Braun et al (61) revealed that germline polymorphisms at AGAP11 were associated with stage-adjusted survival in ovarian cancer. In the current network, hsa-miR-155 was highly expressed, while AGAP11 was expressed at low levels, and may be a prognostic marker in BC. Additionally, 74 DEmRNAs were identified in the current ceRNA network. GO analysis revealed that the main enriched GO biological processes were ‘positive regulation of GTPase activity’, ‘positive regulation of ERK1 and ERK2 cascade’ and ‘transcription from RNA polymerase II promoter’. Analysis of the enriched GO networks of the KM-significant DEmRNAs revealed that nine DEmRNAs were mainly enriched in ‘cytosol’ and ‘positive regulation of epithelial cell proliferation’. Among the nine DEmRNAs, upregulated RASGRP1 expression was demonstrated to be associated with favorable outcome in BC. As a guanine nucleotide exchange factor, RASGRP1 can sensitize Ras and ERK/MAPK cascades. Alterations in RASGRP1 expression have been associated with a variety of diseases, including systemic lupus erythematosus (62), hepatotoxicity (63), and cancer (64,65). Upregulated RASGRP1 expression correlates with improved clinical outcomes in patients with colorectal cancer, and may induce a negative feedback loop to suppress proliferative epidermal growth factor receptor/son of sevenless homolog 1/Ras signaling (64). RASGRP1 may be a novel therapeutic target for the treatment of hepatocellular carcinoma (66). Wang et al (67) reported a positive correlation between RASGRP1 overexpression and survival in triple-negative breast cancer. Some studies have reported that RASGRP1 is associated with the acquired drug resistance of tumors (68,69). In addition, seven other genes identified in the present study (KPNA2, NTRK2, PTGS2, RACGAP1, SFRP1, SHCBP1 ZNF367) have been reported to be associated with the prognosis of tumors (70–75).

Pathway analysis suggested that DEGs were mainly involved in the ‘MAPK signaling pathway’ and the ‘Ras signaling pathway’ in the present study. Additionally, the present results demonstrated that RASGRP1 was enriched in the ‘MAPK signaling pathway’ and the ‘Ras signaling pathway’. MAPKs are a class of serine/threonine protein kinases in cells, and the MAPK signaling pathway is present in most types of cells. MAPK signaling serves an important role in transducing extracellular signals into cells and their nuclei, and induces various cellular biological processes, including cell differentiation, proliferation and apoptosis. Numerous studies have demonstrated that the MAPK signaling pathway is associated with the development of BC (76–79). Yang et al (78) indicated that C-C motif chemokine ligand 28 has been associated with the proliferation and metastasis of BC via the MAPK-mediated anti-apoptotic and metastasis signaling pathways. The androgen receptor pathway is a potential target for the treatment of BC, and has been reported to exhibit crosstalk with the MAPK pathway (77). Previously, the MAPK signaling pathway was demonstrated to be involved in the development of drug resistance in tumors (78). Jia et al (79) reported that Kruppel-like factor 4 (KLF4) could enhance the response of BC cells to tamoxifen by inhibiting the MAPK signaling pathway, indicating that targeting KLF4/MAPK signaling may be a potential therapeutic strategy for treating BC. RAS is one of the most commonly mutated genes linked to human cancer. Anticancer inhibitors targeting the RAS signaling pathway are potential therapeutic agents for the treatment of cancers; however, further investigation is required. The abnormal activation of the Ras/ERK pathway mediates the occurrence and invasiveness of BC, which is notably common in BC. Ras activation is a crucial determinant of the transmission and poor prognosis of ERα+/luminal BC (80,81).

To the best of our knowledge, the present study is the first to report the aberrant expression of INC AC112721.1, AL356479.1, LINC00466 and MIR7-3HG in BC, indicating their potential prognostic value in BC. The present study proposed that LINC00466 and MIR7-3HG interact with hsa-miR-204, which may be involved in the development of BC; to the best of our knowledge, this has not been reported previously. Of note, further bioinformatics analysis of the ceRNA network in BC may aid future investigations. The current results determined that the genes associated with survival in BC were enriched in the ‘MAPK signaling pathway’, which has been reported to be associated with drug resistance. Our group aims to further investigate the pathology of BC via in vivo and in vitro experiments in the future.

Acknowledgements

Not applicable.

Funding

This work was supported by the Science and Education for Health Foundation of Suzhou for Youth (grant nos. kjxw2018030 and kjxw2018032), the Science and Technology Project Foundation of Suzhou (grant no. SS201651), the Education Research Project Foundation of Nanjing Medical University (grant no. FZS-ZD-201701) and the Key Medical Disciplines in Jiangsu Province (grant no. ZDXKC2016007).

Availability of data and materials

The datasets used and/or analyzed during the present study are available from the corresponding author on reasonable request.

Authors' contributions

JJW, YQH and WS carried out data acquisition and processing, and were the major contributors in writing the manuscript. YFL and HW made contributions to data analysis, interpretation and figure editing. WJW and MH directed the whole study, and revised the manuscript. All authors read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References