Construction and Analysis of a Long Non-Coding RNA-Associated Competing Endogenous RNA Network Identified Potential Prognostic Biomarkers in Luminal Breast Cancer

Introduction

Breast cancer (BC) is a major concern, both morbidity and mortality, for women worldwide.¹ Four subtypes of breast cancer (BC) defined by the estrogen receptor (ER), progesterone receptor (PR), and human epidermal growth factor receptor 2 (HER-2) exhibit differences in biological behavior. The molecular mechanism behind them remains unclear. The current target-therapies are mostly concentrated on HER-2 overexpressed subtype, such as Trastuzumab and Pertuzumab,^2,3 while immunotherapies are aimed at triple-negative breast cancer (TNBC) due to its higher immunogenicity.^4,5 Quite a number of prognostic markets and targets were identified.^6–8 Unlike the two subtypes, patients with ER- or PR-positive BC, termed luminal subtype, are considered at a lower risk and endocrine therapy has achieved a considerable success. However, the outcome of luminal BC rather varies after treatment⁹ and drug resistance is urgent to be responded. mTOR inhibitors and CDK4/6 inhibitors showed potential to solve the problem. Their clinical effective suggests that molecular target therapy could also be a wise strategy for luminal BC. Therefore, the identification of prognostic markers, which can be potential targets, is of importance. Recently, genes or proteins were proved to be an oncogene role or promoting metastasis in the luminal subtype,^10–13 while little were verified to be related to overall survival (OS). Moreover, most of them were genomic, along with some microRNAs (miRNA), while Long non-coding RNAs (lncRNA) were rarely mentioned.

LncRNA, a type of RNA with limited protein-coding ability, is previously viewed as transcriptional “noise” without biological functions.¹⁴ Recently, lncRNAs are increasingly reported to be closely involved in progression and metastasis of various cancers.^15–18 In luminal BC, the oncogenic roles of lncRNAs were also unveiled^19,20 but rare, especially specific on luminal BC. And crosstalks between lncRNAs and miRNAs were also identified. For example, lncRNA BLACAT1 was found to promote MCF7 (luminal BC cell) proliferation and metastasis by miR-150-5p/CCR2.²¹ In addition, prognostic roles and relationship with tamoxifen resistance were also mentioned.^13,22

Recently, competitive endogenous RNA (ceRNA), as a whole view of transcription, has gradually become a useful tool for identifying prognostic lncRNAs. The ceRNA hypothesis suggests that lncRNAs can competitively bind to miRNAs, thereby preventing their binding to mRNA, and consequently regulating mRNA expression.²³ The competitive relationship influences the inhibition of miRNAs to mRNAs, resulting in downstream pathological processes like osteoarthritis, cerebral ischemia, and cancers.^24–27 However, the prognosis function of lncRNAs in breast cancer, especially luminal, is less known.

In the study, a ceRNA network was constructed by integrating analysis of several databases and six lncRNAs PART1, IGF2.AS, WT1.AS, OIP5.AS1, SLC25A5.AS1, and AL035706 were identified and verified by human breast cancer specimens. Further, a prognostic model consisted of them was constructed.

Materials and Methods

Download of Data and Pre-Processing

RNA-seq data and mRNA-seq data of BC were obtained from The Cancer Genome Atlas (TCGA) database.²⁸ Samples meet requirements were screened to obtained three matrix files of mRNA, lncRNA, and miRNA expression profiles: 1) female samples; 2) positive for at least one of the Estrogen receptor (ER) or Progesterone receptor (PR) in the cancer samples; 3) both RNA-seq and miRNA-seq data. 4) miRNAs with a ratio of zero values ≤75%.

Relationship with Prediction

The miRcode²⁹ is a database based on the human complete transcriptome of GENCODE annotations to predict miRNA targets, through which we obtained the lncRNA–miRNA relationship files. The miRDB³⁰ is a miRNA target gene prediction database based on high-throughput sequencing experiments. miRTarBase³¹ is a database of integrated, experimentally validated miRNA targets. TargetScanS³² predicts miRNA target genes by searching for conserved 8 mer and 7 mer positions that match miRNA seed sequences point. The three databases are currently used for miRNA target gene predictions. Through this approach, the miRNA–mRNA relationship pair files were screened.

Relationship Screening

Pearson correlation coefficient was calculated between the lncRNA and miRNA in the lncRNA–miRNA relationships, and the miRNA and mRNA in the miRNA–mRNA relationships. The significant negative correlation (p<0.05) relationship pairs, which served as a pre-selected ceRNA regulatory relationship, were screened.

Building a ceRNA Network

Files that were only in the lncRNA–miRNA or miRNA–mRNA relationship files were excluded and the remained lncRNA–miRNA and miRNA–mRNA files were merged, and visualized by the Cytoscape software.

mRNA Enrichment Analysis in the ceRNA Network

KEGG is a group of databases ranging from genomes, biological pathways, and associated diseases, even drugs and chemical materials.³³ Gene ontology analysis (GO) is commonly used for interaction network analysis of genes and their RNA or protein products.³⁴ By KEGG and GO, mRNA enrichment analysis was performed in the ceRNA network.

Topological Analysis and Stability Analysis of the ceRNA Network

The topology properties of the network were analyzed using the NetWorkAnalyzer toolKit 10 in Cytoscape9.^35,36 The NetWork Analyzer mainly includes network distribution, number of connections, average aggregation coefficient, etc. In the present study, we calculated the number of connections, the length of the path, and the proximity center of the node.

Prognostic Analysis of ceRNA Modules

The clinical information of the above samples was download from the TCGA database, and then the survival data of each sample were extracted, which combined the previously obtained expression spectrum data to perform a Cox survival analysis on the lncRNA nodes in the network then identify prognostic lncRNAs.

Single-Factor Regression, Multivariate Cox Regression Analysis, and Survival Analysis

Multivariate Cox regression analysis and survival analysis were performed on prognostic lncRNAs by SPSS 19.0, and then a prognostic model was constructed and survival curves were plotted.

Patients and Tissue Samples

Paired tumor specimens and their adjacent non-tumor tissues were derived from 20 luminal A patients who underwent surgical resection at Second Affiliated Hospital, Zhejiang University School of Medicine (Hangzhou, China) between Jan 2019 and May 2019. All patients recruited in the study received no pre-operative chemotherapy or radiotherapy. Specimens were frozen in liquid nitrogen immediately and stored at −80°C. This study was approved by the Ethics Committees of Second Affiliated Hospital, Zhejiang University School of Medicine.

RNA Extraction and RT-qPCR Analysis

Using Trizol Reagent (Invitrogen), total RNA was extracted from tissues. A total of 2 μg RNA was reverse transcribed to cDNA according to the PrimeScript RT reagent Kit (TaKaRa, Dalian, China). Then, qPCR assay was performed by the SYBR Premix Ex Taq™ II kit (TaKaRa), primers, and cDNA templates on the ABI PRISM 7900 Sequence Detection System (Applied Biosystems, Foster City, CA, USA). The 2^−ΔΔCt method was used for data quantification. GAPDH was set as the internal control. The sequences of six lncRNAs were as follows

IGF2-AS: Forward primer CCTGCCTAGAGCTCCCTCTT;

Reversed Primer CCATCCTCACCCAGGAACAG.

WT1-AS: Forward primer CTTTCGACTAGCGCCTCTCC;

Reversed primer CTTTCGACTAGCGCCTCTCC.

OIP5-AS1: Forward primer CCTCATGCAGTGCCATCTGA;

Reverse primer TATCCACCTTGGGTTGCAGG.

SLC25A5-AS1: Forward primer AGCGGCATCTGTCATGTTGA;

Reverse primer AAATCGGCCATTTGCTTCGC.

AL035706.1-201: Forward primer GCCAAAGTGCTGGGATTACA;

Reverse primer GCCACGGGAATCAAATAAGA.

PART1: Forward primer AAACGCCTGAGGACTGAGAA

Reverse primer TCTCCTGCTTGCCAAATCTT

Logistic Regression Analysis for Pathological Parameters

Using Ki67 (1 ≤ 10%, 0 = ≤10%), Grade (1 = Ⅲ, 0 = Ⅰ or Ⅱ), diameter (1 ≤ 1.5, 0 = ≤1.5), age (1 ≤ 50, 0 = ≤50), pathological types (1 = invasive, 0 = non-invasive), and the surgical strategy (1 = conservation, 0 = radical) as the dependent variables and the expression of six prognostic lncRNAs as the independent variable to perform binary logistic regression analysis. P < 0.05 is considered as a statistically significant difference. Stata 12.0 was used for statistical analysis.

Results

lncRNA–miRNA Relationship

From TCGA, the RNA-seq data contained 113 normal samples and 1102 cancer samples and the miRNA-seq data contained 104 normal samples and 1096 cancer samples were download. After screening, we got 884 samples, including 102 normal samples and 782 cancer samples. Subsequently, we extracted the expression profiles of mRNA and lncRNA from the matrix files obtained from the RNA-seq data. An expression profile of 14,166 lncRNAs was revealed, and by miRcode, we further identified 7098 pairs of lncRNA–miRNA, including 1449 lncRNAs and 35 miRNAs.

miRNA–mRNA Relationship

An expression profile containing 2253 miRNAs was collected. After excluding miRNAs with a ratio of zero values >75% in all samples, 875 miRNAs remained. Using miRDB, miRTarBase, and TargetScan to predict the target genes of these miRNAs, we obtained 19,616 pairs, including 5322 target genes simultaneously predicted by all three databases and 785 miRNAs.

Relationship Filtering

By screening, we obtained 295 pairs of lncRNA–miRNA and 3204 pairs of miRNA–mRNA. After a future filtering of the non-compliant relationship, we finally identified 240 pairs of lncRNA–miRNA and 115 pairs of mRNA–miRNA. Several of those were shown in Tables 1 and 2.

Table 1 Several of the lncRNA–miRNA Pairs Obtained from This Analysis

Table 2 Several of the miRNA–mRNA Pairs Obtained from This Analysis

Construction and Visualization of the ceRNA Network

After combining the aforementioned 240 pairs of lncRNA–miRNA and 115 pairs of mRNA–miRNA, we used the Cytoscape software for network visualization to obtain a ceRNA network map. The network map contained 356 edges and 288 nodes, including 160 lncRNAs, 15 miRNAs, and 113 mRNAs (Figure 1).

Figure 1 Diagram of the ceRNA network. Red point: miRNA; blue point: lncRNA; yellow point: mRNA. The size of the node point represents the degree of the node. Higher degrees indicate larger points.

Topological Analysis and Stability of the ceRNA Network

Topological analysis showed the distribution of the degree of nodes (Figure 2A). The degree of most nodes was low. Only a small number showed a high degree, suggesting that these nodes acted as a hub for the entire network.

Figure 2 Topological analysis and stability analysis of the ceRNA network. (A) Node degree distribution map. (B) Closeness Centrality distribution map. (C) Shortest path distribution map. (D) The density map of node degree distribution.

The Closeness centrality (CC) of the node could be used to calculate the number of connection steps between nodes. The more concentrated nodes yield higher scores; thus, the CC indicated the shortest path. Figure 2B illustrated that nodes with a relatively low number of connections were relatively concentrated, whereas those with higher connectivity were scattered.

The path reflected the combination of all nodes in the network. Figure 2C showed the shortest path length distribution, suggesting that the path length distribution was concentrated, as the extreme values at both ends were lower. The upper and lower limit was 8 and 1, respectively, indicating that most nodes in the network could be connected through a shorter path.

The density map of node degree distribution is shown in Figure 2D. We found that node density decreased sharply with an increase in node degree, suggesting that most nodes in the network were isolated. In the process of disease, only a few key nodes changed and interacted with adjacent nodes. Subsequently, co-expression occurred, which in turn triggered the downstream biological processes.

Enrichment Analysis

In order to observe the function of the constructed ceRNA network, we selected all the mRNAs included in the network to perform functional enrichment analysis.³⁷

The results of the Gene Ontology (GO) analysis of the ceRNA network showed that most of the enriched biological processes were related to blood agglutination and hormones, suggesting that hormone regulation may be the main processes, by which key RNAs involved in the development of BC (Figure 3A).

Figure 3 Results of the enrichment analysis. (A) GO and (B) KEGG.

The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis revealed that choline metabolism, T-cell receptor signaling, the tumor necrosis factor signaling pathway, and the mitogen-activated protein kinase signaling pathway were related to the key nodes of the ceRNA (Figure 3B).

lncRNA Node Survival Analysis

The 160 lncRNAs included in the network were subjected to Cox survival analysis, and six lncRNAs significantly correlated with OS were identified (P<0.05); PART1, IGF2.AS, WT1.AS, OIP5.AS1, and SLC25A5.AS1 were associated with an adverse OS (HR>1) while AL035706.1 with a good one (HR<1) (Table 3).

Table 3 Prognosis-Related lncRNAs in the ceRNA Network

Prognosis-Related lncRNA Molecules as Potential Markers

We performed single-factor regression on the six identified lncRNAs, and the forest map was shown in Figure 4. The high expression of five of those was associated with poor prognosis, whereas the low expression level of the remaining one lncRNA was associated with poor prognosis.

Figure 4 Forest map of the six identified prognosis-related lncRNAs.

Different Expression of Six lncRNAs Between Normal and Tumor Tissue of Breast Cancer Patients

We next verified six lncRNAs expression by clinical samples. Ten pairs of normal and tumor tissue from fresh human breast cancer specimens (ER+) were collected and examined the expression of six lncRNAs. The results showed that among five adverse-OS ones, PART1, IGF2.AS, WT1.AS, and SLC25A5.AS1 significantly higher expressed in tumor than normal, especially for PART1. While OIP5.AS1 expressed no difference between normal and tumor. AL035706.1 was higher expressed in normal, which was consistent with the former results (Figure 5).

Figure 5 The expression of the six identified prognosis-related lncRNAs between normal and tumor tissue of luminal breast cancer specimens collected from our center. N=20; *P<0.05, **P<0.01, ***P<0.001.

The Expression of Six lncRNAs Correlates with Higher Pathological Parameters

To evaluate whether six lncRNAs’ expression correlated with pathological parameters or not, we chose the Ki-67, tumor grades and tumor diameters of patients to help perform the logistic regression analysis. We found that the expression of PART1 was significantly in accordance with all three parameters, and AL035706.1 was adverse with three parameters but had not reached a significant difference. The higher expression of WT1.AS and SLC25A5.AS1 was associated with an advanced grade and the higher expression of OIP5.AS1 with a higher Ki-67. However, IGF2.AS did not show a relation with the parameters in our samples. To exclude the other clinical parameters’ influence, we also evaluated the age, pathological types, and surgical strategies, and the results were negative (Table 4).

Table 4 The Correlation Between Prognosis-Related lncRNAs and Pathological Parameters in Clinical Samples

Prognostic Model of lncRNA

Further, we used the expression profiles of these six lncRNAs in a multivariate Cox regression analysis to construct a prognostic model: the Risk Score = (0.0906×PART1) + (0.1506×IGF2.AS) + (0.2048×WT1.AS) + (0.2087×OIP5.AS1) − (0.8068×AL035706.1) + (0.2072×SLC25A5.AS1). Subsequently, we calculated the risk coefficient by taking the expression profile of each lncRNA into the model. The samples were classified according to the median risk coefficient. As shown in Figure 6, the high-risk group had significantly more deaths than the low-risk group. Moreover, a higher risk coefficient was linked to a higher expression level for lncRNAsPART1, IGF2.AS, WT1.AS, OIP5.AS1, and SLC25A5.AS1. The opposite was observed for AL035706.1.

Figure 6 The relationship between the expression of the six identified lncRNAs and risk scores. The horizontal axis represents the risk coefficient score of the sample, increasing from left to right.

As shown in Figure 7A, we performed a survival analysis and plotted Kaplan–Meier curves to determine prognostic differences between the high- and low-risk samples of the model. A significant difference in prognosis was observed between the two groups. In order to verify the heterogeneity and stability of the model, we analyzed the Receiver operating characteristic (ROC) curve of the model, as shown in Figure 7B. The prognostic model exhibited a markedly greater area under the curve, and prognosis was classified as high and low risk according to these six lncRNAs. These findings suggested that the six identified lncRNAs may serve as meaningful prognostic markers.

Figure 7 (A) Kaplan–Meier curves from the multivariate regression analysis for the prognostic classification of the six lncRNAs. (B) ROC curve from the multivariate regression analysis for the six lncRNAs.

Discussion

Luminal BC is the most common type of whole BC.⁹ Owing to advanced screening and endocrine therapy, some success in the treatment of luminal BC has recently been achieved. However, the mortality rate has not declined, partly because of a lack of efficient biomarkers. Considering the existed heterogeneity of luminal BC reflected by its inconsistent clinical outcomes, it is necessary to discover effective biomarkers, especially for luminal BC. Previous researches have shown several dysregulated genes that involved in the progression of BC and possessed great potential to be used as biomarkers. However, these markers focused on TNBC or distant metastasis.^{6–8,10–13} Recent studies have revealed that lncRNAs played vital roles in the development of cancers.^21,22 Further, germline mutations can also predict tumor recurrence in breast cancer patients.³⁸ In the present study, we used ceRNA to detect and validate potential prognostic lncRNAs in the luminal BC.

Firstly, using TCGA database, we obtained the RNA-seq data and miRNA-seq data of BC. These data were compared with the data in the miRcode²⁹ database to get 7089 pair files of the lncRNA–miRNA relationship. Subsequently, we defined the miRNA-targeted genes that were predicted by all three databases (miRDB, miRTarBase and TargetScan) as the final target genes, thereby obtaining 19,616 pair files of the miRNA–mRNA relationship. After screening, 240 pairs of lncRNA–miRNA and 115 pairs of mRNA–miRNA remained. Thus, a ceRNA topology network with 356 edges and 288 nodes, including 160 lncRNAs, 15 miRNAs, and 113 mRNAs was constructed. Through topological analysis, we observed that only a small number exhibited a high degree, suggesting that these act as hubs for the entire network. Most nodes in the network can be connected through shorter path correlation and showed isolation. This result indicates that only a few co-expressed key nodes change and interact with adjacent nodes in the process of BC, which in turn triggered the downstream biological processes. Subsequently, we used the GO and KEGG databases to perform enrichment analyses to elucidate involved biological process and pathways.³⁷ The analyses suggested that hormone regulation, choline metabolism, T-cell receptor signaling, the tumor necrosis factor (TNF), and the mitogen-activated protein kinase (MAPK) signaling pathway were involved in the development of BC. After Cox survival analysis, six lncRNAs, namely PART1, IGF2.AS, WT1.AS, OIP5.AS1, AL035706.1, and SLC25A5.AS1, significantly correlated with the prognosis were finally identified. According to single-factor regression and forest map analysis, the high expression level of PART1, IGF2.AS, WT1.AS, OIP5.AS1, and SLC25A5.AS1 was associated with poor prognosis, while the low expression of AL035706.1 was linked to poor prognosis. Furthermore, we used the expression profiles of these six lncRNAs in a multivariate Cox regression analysis to construct a prognostic model. The survival analysis and Kaplan–Meier curves showed a significant difference in prognosis between the high- and low-risk groups classified by the model, and a markedly greater area under the curve verified the heterogeneity and stability of the model.

More importantly, 20 normal and tumor paired luminal breast cancer specimens were collected to examined six prognostic lncRNAs’ expression and their correlation with clinical pathological parameters. We verified that five poor-prognostic lncRNAs were higher expressed in tumor than normal, and their expression, except for IGF2.AS, was associated with, to a degree, adverse parameters, such as Ki-67, grade, and tumor diameter. The results supported that these lncRNAs had a specification in tumor so that could be candidate prognostic markers. As for the positive-prognostic one, AL035706.1 was higher expressed in the normal tissue and showed a tendency, but not a significance, to oppose adverse parameters. For this, we supposed a larger number of specimens could give a definite conclusion.

We found six potential lncRNAs, namely PART1, IGF2.AS, WT1.AS, OIP5.AS1, AL035706.1, and SLC25A5.AS1. Only the high expression of AL035706.1 predicted a better prognosis, others were all with worse. Consistent with our findings, it was previously reported that PART1 could predict early recurrence in patients with hepatocellular carcinoma after curative resection, and modulate the toll-like receptor (TLR) pathway to influence cell proliferation and apoptosis in prostate cancer cells.^39–41 The research showed that IGF2.AS was related to metabolism, and inhibition of IGF2.AS promoted angiogenesis in patients with type 2 diabetes.⁴² While the other four lncRNAs have not been well investigated thus far and all of their prediction in OS were not mentioned. We proposed that they may be promising new targets for further researches.

Conclusion

In the present study, we comprehensively analyzed several databases to construct a ceRNA network in luminal BC and identify six OS-related lncRNAs, by which a prognosis prediction module was constructed. The six lncRNAs’ higher expression and pathological features were verified by our paired breast cancer specimen. Therefore, we proposed that these six lncRNAs can be useful biomarkers for predicting the prognosis of patients with luminal BC.