Construction and Analysis of the Dysregulated ceRNA Network and Identification of Risk Long Noncoding RNAs in Breast Cancer

Breast cancer (BRCA) is the second leading cause of cancer-related mortality in women worldwide. However, the molecular mechanism involved in the development of BRCA is not fully understood. In this study, based on the miRNA-mediated long non-coding RNA (lncRNA)–protein coding gene (PCG) relationship and lncRNA–PCG co-expression information, we constructed and analyzed a specific dysregulated lncRNA–PCG co-expression network in BRCA. Then, we performed the random walk with restart (RWR) method to prioritize BRCA-related lncRNAs through comparing their RWR score and significance. As a result, we identified 30 risk lncRNAs for BRCA, which can distinguish normal and tumor samples. Moreover, through gene ontology and Kyoto Encyclopedia of Genes and Genomes pathway analysis, we found that these risk lncRNAs mainly synergistically exerted functions related to cell cycle and DNA separation and replication. At last, we developed a four-lncRNA prognostic signature (including AP000851.1, LINC01977, MAFG-DT, SIAH2-AS1) and assessed the survival accuracy of the signature by performing time-dependent receiver operating characteristic (ROC) analysis. The areas under the ROC curve for 1, 3, 5, and 10 years of survival prediction were 0.68, 0.61, 0.62, and 0.63, respectively. The multivariable Cox regression results verified that the four-lncRNA signature could be used as an independent prognostic biomarker in BRCA. In summary, these results have important reference value for the study of diagnosis, treatment, and prognosis evaluation of BRCA.

Introduction

Breast cancer (BRCA) is one of the most prevalent malignancies and is the second leading cause of cancer-related mortality in women worldwide (Bray et al., 2018; DeSantis et al., 2019). Although there have been several advancements in both surgical and adjuvant therapy, BRCA remains a significant threat to female health due to high incidence and poor prognosis. There is an urgent need for novel and effective biomarkers for clarifying mechanism of early BRCA and providing therapeutic targets for BRCA patients (Chi et al., 2019).

Long non-coding RNAs (lncRNAs) are a group of RNAs with length >200 bp, which serve as key regulators in diverse cellular functions such as development, differentiation, and apoptosis (Ulitsky and Bartel, 2013; Quinn and Chang, 2016; Trovero et al., 2020). The important function of lncRNA is that it can act as a competitive endogenous RNA (ceRNA) to regulate the expression level of other transcripts especially for protein coding gene (PCG) by sponging miRNA (Quinn and Chang, 2016). Growing evidences demonstrated that lncRNAs had been indicated as important molecules in tumorigenesis (Gibb et al., 2011). Herrera-Solorio et al. (2020) showed that lncRNA SOX2-OT modulates an orchestrated resistance mechanism, promoting poor prognosis and human lung malignancy through genetic, epigenetic, and posttranslational mechanisms. Xiao et al. (2020) found that the lncRNA MAFG-AS1, which is highly expressed in bladder urothelial carcinoma, is correlated with aggressive characteristics and poor prognosis of bladder urothelial carcinoma. As for BRCA, IRNAS HOTAIR, SPRY4-IT1, GAS5, MATAR25, PANDAR, and MATAR25, new players in tumor development and prognosis, have shown owning potential clinical applications in BRCA (Soudyab et al., 2016; Nagini, 2017; Chang et al., 2020). In addition, viable ways have been considered to predict the potential BRCA lncRNAs by performing high-throughput data based on bioinformatics methods (Guo et al., 2018; Chi et al., 2019; Wang et al., 2019).

In our study, by analyzing the BRCA expression profile from The Cancer Genome Atlas (TCGA) and lncRNA-related databases, we constructed and analyzed a specific dysregulated BRCA-associated lncRNA-PCR ceRNA network and performed the random walk with restart (RWR) method to prioritize BRCA-related lncRNAs through comparing their RWR score and significance. At last, we identified 30 risk lncRNAs associated with BRCA and constructed a prognostic signature based on the TCGA expression data with clinical survival characters. This study has important reference value to accelerate the discovery of molecular biomarkers for the study of diagnosis, treatment, and prognosis evaluation of BRCA.

Materials and Methods

Datasets Across Breast Cancer

The PCG and lncRNA expression profiles of BRCA with FPKM values were obtained from the UCSC Cancer Browser¹, which provided an open-access portal to download data from TCGA. In total, we acquired 1,222 BRCA samples that were involved in 1,075 patients with complete clinical follow-up information. We performed a two-step filter for PCG and lncRNA expression profiles; the aim is to ensure detection reliability and reduce noise. First, we extracted only paired patient samples with tumor and adjacent nontumor tissue for differentially analysis. Second, lncRNAs or PCGs with an average expression value of less than 1 were removed in the tumor and adjacent nontumor tissue. Finally, we obtained 224 tumor and adjacent nontumor tissue from 112 patients, including 1,251 lncRNAs and 13,356 PCGs. Processed gene expression data and clinical data were provided in Supplementary Material. For the calculation of Pearson correlation coefficient (PCC), log2 transformation was performed to lncRNA/TF expression profiles with raw expression values.

Breast Cancer-Associated Known PCGs

We downloaded BRCA-related PCGs from DisGeNET (v7.0) (Piñero et al., 2020). DisGeNET is a discovery platform containing one of the largest publicly available collections of genes associated with human diseases, which integrates data from expert curated repositories and the scientific literature (Piñero et al., 2020). We extracted 318 PCGs associated with BRCA from DisGeNET (Supplementary Table 1).

Identification of the Differentially Expressed PCGs/LncRNAs

Fold change and statistical significance were computed for each PCG/lncRNA in expression profiles by limma package, which is a common, effective R/Bioconductor software package for differential expression analyses (Ritchie et al., 2015). The lncRNAs and PCGs with P < 0.05 or | log2 fold change| > 1 were considered to be differentially expressed (DE) lncRNAs/PCGs.

Construction of a Specific LncRNA–PCG ceRNA Network for Breast Cancer

We performed two steps to construct specific LncRNAs-PCGs ceRNA network for breast cancer (SLGCeNBC). Firstly, we calculated the correlation of co-expression between DE PCGs and lncRNAs using the PCC method. PCC can be used to measure the linear relationship between lncRNA and PCG expression (Zhang et al., 2018). In this study, we considered that the lncRNA–PCG pairs with PCC > 0.5 and P < 0.01 showed a potential expression correlation. All lncRNA–PCG pairs meeting the threshold were merged into the lncRNA–PCG co-expression network.

Secondly, we performed the hypergeometric test to identify miRNA-mediated lncRNA–PCG pairs. Previous studies proposed and demonstrated that the number of common targeting miRNAs between lncRNAs and PCGs determined the ceRNA cross-talk strength. More common targeting miRNAs could produce a more strength ceRNA cross-talk pair. Thus, we downloaded the interaction of miRNA–PCG/lncRNA from StarBase (Li et al., 2014), mirTarbase (Chou et al., 2018), TargetScan (Grimson et al., 2007), LncBase (Karagkouni et al., 2020), and MiRcode (Jeggari et al., 2012) databases. We downloaded 714,288 miRNA–PCG interaction pairs and 766,809 miRNA–lncRNA interaction pairs, including 13,295 PCGs, 13257 lncRNAs, and 2,593 miRNAs. Then, we used these data to perform the hypergeometric test (Feng et al., 2019). We considered P-value < 0.01 as statistically significant. The P-value was measured as the following:

$p=1-\sum _{i=0}^{r-1}\frac{\left(\begin{array}{c}\hfill t\hfill \\ \hfill i\hfill \end{array}\right)\left(\begin{array}{c}\hfill m-t\hfill \\ \hfill n-i\hfill \end{array}\right)}{\left(\begin{array}{c}\hfill m\hfill \\ \hfill n\hfill \end{array}\right)}$

where m stands for the total number of human genome miRNAs, t stands for the number of miRNAs interacting with the PCG, n stands for the number of miRNAs interacting with the lncRNAs, and r stands for the number of miRNAs shared between PCGs and lncRNAs.

Random Walking Analysis

Here, we performed RWR to determine ranking for BRCA-related lncRNAs. A random walk in network was defined as an iterative walker’s transition from its certain node to a randomly selected neighbor that started from a given node (e.g., “PCG x” was a known PCG associated with BRCA) (Zhang et al., 2018). The random walk performed had capacity of restart with probability r in every time step at node “PCG x.” The random walk with restart was defined as the following:

$p_{t+1=}\left(1-r\right)Wpt+rp0$

where W represents the column-normalized adjacency matrix of the network, p_t is a vector whose size is equivalent to the number of nodes in the network, and the i-th element holds the probability of being at node i at time step t.

The initial probability vector p₀ was constructed such that 1 was assigned to the nodes representing known PCGs associated with BRCA, and other nodes with 0. We considered that the role of PCGs related to disease was equivalent in the network. Vector p would be in the steady state at time step t, where t approached infinity as a limit. The iteration would be finished till the change between p_t and p_t+1 falls below 10^–10.

We scored for each lncRNA to prioritize lncRNAs related to BRCA by performing the RWR algorithm in a specific LncRNA–PCG ceRNA network for breast cancer (SLGCeNBC) and performed statistical significance analysis for the score of every lncRNA. The statistical significance was determined by comparing the scores of lncRNAs in the network following n iterations with SLGCeNBC perturbation. To maintain the network topological properties, random sampling without replacement was performed when doing the random disturbance. In iterations, the times that the score of every lncRNA was higher than the real one were recorded as m. The P-value for every lncRNA was the ratio of m and n. In this study, n was set at 10,000 times.

Constructing the Prognostic lncRNA Signature

To identify the best score cutoff of the selected lncRNAs for grouping patients most significantly, we employed for optimal cutoff identification by using the R package “maxstat” (Hothorn and Zeileis, 2008). The survival outcomes of the two groups were estimated by Kaplan–Meier analysis (Bland and Altman, 1998). Then, we used the survival-related lncRNAs (log-rank test P < 0.05) above to perform Cox regression analysis (Cox, 1972) and construct risk models as follows:

$\mathrm{Risk}\mathrm{Score}=\sum _{i=1}^n\left(x_i{e}_i\right)$

where n is the number of lncRNAs, e_i is the expression value of the lncRNA_i, and x_i is the coefficient of lncRNAs in Cox regression analysis. Finally, we used the signature with the minimum log-rank P-value as the best prognostic marker (Guo et al., 2016). Cox regression analysis was performed to explore the predictive independence of the lncRNA signature. R software² with R packages including timeROC, survival was used for statistical analysis, where a P-value of < 0.05 was considered statistically significant.

Results

Construction of a Specific LncRNA–PCG ceRNA Network for Breast Cancer

Based on the interaction of miRNAs–PCGs/lncRNAs from the public database, we merged all miRNA–PCG pairs and miRNA–lncRNA pairs and then obtained a global PCG–miRNA–lncRNA triple network. PCG and lncRNA which shared at least one common miRNA were reserved. Then, we performed the hypergeometric test to identify PCG–lncRNA pairs with P-value < 0.01, which yielded more than candidate 800,000 lncRNA–PCG interactions. All the pairs were merged into a miRNA-mediated lncRNA–PCG network.

Second, based on the expression profile of 112 pairs of tumor and adjacent nontumor tissue, we obtained DE lncRNAs/PCGs of BRCA by limma with P < 0.05 and | log2 fold change| > 1, which produced 316 and 2,463 DE lncRNAs and PCGs, respectively. Then, we calculated the Pearson correlation between DE PCGs and lncRNAs with PCC > 0.5 and P < 0.01. All the co-expressed lncRNA–PCG pairs were merged into a co-expression network. As for the lncRNA–PCG co-expression network, it contained 307 lncRNAs, 2,105 PCGs, and 57,216 co-expression relationships. Obviously, it is specific for BRCA.

Finally, we extracted common lncRNA–PCG pairs in the above two networks to construct SLGCeNBC. We defined the PCGs or lncRNAs as the nodes of the network. The whole identification process is shown in Figure 1. In total, 259 lncRNAs, 1,384 PCGs, and 21,702 edges were included in SLGCeNBC (Figure 2A and Supplementary Table S2).

FIGURE 1

Figure 1. Schematic overview in this study. (i) LncRNA–PCG ceRNA network. PCG and lncRNA shared at least one common miRNA and considered that a P-value < 0.01 was significant by applying the hypergeometric test. (ii) LncRNA–PCG co-expression network. The correlation of co-expression between DE PCGs and lncRNAs using the PCC (PCC > 0.5 and P < 0.01). (iii) Extracted common lncRNA–PCG pairs in (i) and (ii) to structure SLGCeNBC. (iv) The BRCA-associated known PCGs (seed nodes) were mapped into the SLGCeNBC, and the RWR method was performed on this network. Finally, we the ranked candidate lncRNAs according to the steady probability of RWR.

FIGURE 2

Figure 2. SLGCeNBC and three lncRNAs appeared in two dimensions. (A) SLGCeNBC. The red and blue nodes represent lncRNAs and PCGs, respectively. A lncRNA and PCG were connected by an edge if they had the ceRNA and co-expression relationships (Supplementary Table 2). (B) The true node degree distribution of SLGCeNBC; the degree distribution of all nodes followed the power law distribution approximately with a slope of –1.054 and R² = 0.839. (C) The lncRNAs (MAGI2-AS3, MIR99AHG, and LINC00641) are ranked in the top 20 in at least two dimensions (degree, betweenness, and closeness). (D) GO enrichment analysis for lncRNAs MAGI2-AS3, MIR99AHG, and LINC00641 through their related PCGs that share at least two of three lncRNAs in SLGCeNBC. The X-axis and bubble color are the –log10 of P-value; the Y-axis is the names of the GO terms. Bubble size indicates the number of PCGs annotated to the GO term.

Topological Analysis of SLGCeNBC

The degree distribution of all nodes followed the power law distribution approximately with R² = 0.839 (Figure 2B). A handful of nodes with a high degree in the networks were defined as hubs that linked many nodes; most nodes in networks had few interactions. We firstly analyzed the topological properties of the SLGCeNBC and calculated the degree, closeness, and betweenness of the SLGCeNBC, respectively. We ranked all the nodes’ topological features of the network and listed the top 20 nodes of degree, betweenness, and closeness, respectively (Table 1). Interestingly, we found that three lncRNAs (MAGI2-AS3, MIR99AHG, and LINC00641) appeared in at least two dimensions (Figure 2C). For the lncRNA MAGI2-AS3, we found that it significantly downregulated (log2FC = −2, P = 1.74E-38) in differential expression analysis. Some studies showed that overexpression of MAGI2-AS3 in BRCA cells MCF-7 would inhibit the Wnt/β-catenin pathway and inhibit cell proliferation and migration. MAGI2-AS3 may act as a cis-acting regulatory element downregulating the DNA methylation level of the MAGI2 promoter region (Du et al., 2019; Xu et al., 2021). Meng et al. (2020) verified that overexpression of MIR99AHG promoted gastric cancer cell proliferation and invasion via the miR577/FOXP1 axis. Other experiments have shown that MiR-577 inhibits EMT and metastasis of BRCA by targeting RAB25 (Yin et al., 2018). From these results, it suggested that MIR99AHG may have an effect on BRCA via the miR577/RAB25 axis, which provided a suggestion for further experiments. For LINC00641, Mao et al. (2020) confirmed that LINC00641 inhibits BRCA cell proliferation, migration, and invasion by sponging miR-194-5p.

TABLE 1

Table 1. The top 20 lncRNAs/PCGs in degree, betweenness, and closeness.

To further explore the function of the above three lncRNAs, we performed Gene Ontology (GO) enrichment analysis for the lncRNAs through their related PCGs that shares at least two of three lncRNAs in SLGCeNBC. The result of the GO biological process contained “positive regulation of endothelial cell proliferation,” “cell adhesion,” “angiogenesis,” “regulation of small GTPase mediated signal transduction,” “positive regulation of endothelial cell migration,” “axon guidance,” “phosphatidylinositol phosphorylation,” “positive regulation of apoptotic process,” “regulation of phosphatidylinositol 3-kinase (PI3K) signaling,” and “steroid hormone mediated signaling pathway” (Figure 2D). Many researches had shown that these biological processes were closely associated with BRCA. For example, some deregulation Arf isoforms from the small GTPase subfamily induce cancer formation and progression by enhancing cell proliferation through the activation of mitogen-activated protein kinases (MAPK) and ribosomal protein S6 kinase beta-1 (p70S6K) (Davis et al., 2016; Li et al., 2017). High-level amplification of ARF1 from the Arf subfamily is associated with increased PCG expression and poor outcomes of patients with BRCA (Xie et al., 2016). Overexpression of Ras from the small GTPase subfamily has been found in more than 15% of human tumors (Goitre et al., 2014). The study indicated that Ras which is upregulated in BRCA can promote BRCA cell proliferation, migration, and invasion due to their capability to alter integrin-mediated cell adhesion (Di et al., 2015). For the biological processes, PI3K is the most common altered pathway in ER-positive BRCA and PI3K/AKT is one of the most critical signal pathways for cancer (Hamilton and Infante, 2016; Devanathan et al., 2020). These results recommended that the three lncRNAs with higher degree, betweenness, and closeness were important in the network and played a crucial role in the origin and development of BRCA.

Identifying Risk lncRNAs by Random Walk With Restart

We mapped 318 PCGs associated with BRCA from DisGeNET (v7.0) (Piñero et al., 2020) into SLGCeNBC. The result showed that there are 40 PCGs mapped into SLGCeNBC. The 40 PCGs acted as the seed nodes (Supplementary Tables 1, 3), and the method of RWR (see section “Materials and Methods”) was performed to prioritize BRCA highly related lncRNAs. The initial score of the seed nodes was set at 1, and the scores of all lncRNA node were calculated. To establish whether the lncRNA scores were significantly higher than the random case, we perturbed SLGCeNBC and performed the RWR 10,000 times. As a result, we identified 30 lncRNAs whose scores were significantly higher than those of the random case (P < 0.05, Table 2 and Supplementary Table 4). Here, because the significant results were produced by inputting known BRCA genes and 10,000 times network permutations, these genes were located in the neighbors of the disease genes and considered as the potential synergetic regulators of the disease genes. Thus, all 30 lncRNAs were considered to be risk lncRNAs of BRCA. We showed that the real scores for risk lncRNAs from RWR were higher than the scores for the non-risk lncRNAs (P = 2.266e-06, Wilcoxon rank-sum test). Thirteen of 30 risk lncRNAs have been reported to be associated with the occurrence, progression, and survival of tumor. Particularly, some studies had shown that the nine risk lncRNAs were closely related to the proliferation, metastasis, and survival of BRCA (Table 2). For instance, Zhang et al. (2017) revealed that AP000439.3 could regulate the expression of CCND1 through enhancing estrogen receptor induction of CCND1 and function as a key regulator of the cell cycle in BRCA. Ji et al. (2020) demonstrated that LINC00665 promoted BRCA progression and induced an epithelial–mesenchymal transition-like phenotype via the upregulation of LIN28B expression (Ding et al., 2020).

TABLE 2

Table 2. Risk lncRNA information.

We performed bidirectional hierarchical clustering to further investigate the risk lncRNAs. In the heatmap (Figure 3A), we found that the lncRNAs classified the samples into adjacent nontumor tissue and tumor tissue, suggesting that these lncRNAs possessed potential for diagnosis and therapy of BRCA. Thereafter, we counted the amount of risk lncRNAs that interacted with each PCG. We found that the top two PCGs were ESR1 and PARD6B, which connected with 11 and 10 risk lncRNAs, respectively (Figure 3B). The evidence indicated estrogen receptor-alpha (ERalpha) binding to all identified SRC-3 genomic binding sites from E2-treated cells and confirmed the ability of SRC family coactivators to regulate the expression of one of these PCGs, PARD6B/Par6 (Labhart et al., 2005). We also counted the number of the PCGs that interacted with risk lncRNAs. We found that the top two risk lncRNAs were TYMSOS and AC092718.4, each of which connected with 103 PCGs (Figure 3C). We performed GO enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis for these PCGs related to two risk lncRNAs, respectively. We found that TYMSOS and AC092718.4 were mainly enriched in the biological processes and pathways related to the cell cycle, such as “cell division” and “proliferation.” Interestingly, the result of TYMSOS was similar as that of AC092718.4 (Figure 3D and Supplementary Table 5, Jaccard similarity coefficient = 0.52). However, the similarity of PCGs which was regulated by TYMSOS and AC092718.4 was significantly lower than that of enrichment results (Figure 3C, Jaccard similarity coefficient = 0.35). These results suggest that risk lncRNAs may regulate in coordination the occurrence and development of BRCA.

FIGURE 3

Figure 3. Analysis of the identified risk lncRNAs. (A) The heatmap of 30 risk lncRNAs based on their expression. The columns represented 224 samples of 112 BRCA patients, and the rows represented lncRNAs. The risk lncRNAs classified the samples into tumor tissue and adjacent nontumor tissue. (B) Top two PCGs with high lncRNA interactive relationships. Red nodes represent lncRNAs, and blue nodes represent PCGs. (C) Top two risk lncRNAs with high PCG interactive relationships. (D) GO enrichment analysis for lncRNAs TYMSOS and AC092718.4 through their related PCGs, respectively. The X-axis shows the –log10 of the P-value; the Y-axis shows the names of the GO term. The number at the top of the bar indicates the number of PCGs annotated to the GO term.

Synergistic Regulation of Risk lncRNAs

To further investigate the synergistic regulation function for risk lncRNAs, we focused on the PCGs regulated by more than three risk lncRNAs (Figure 4A) and performed GO function and KEGG pathway enrichment analysis. The results showed that the PCGs were significantly enriched to the KEGG pathway containing “Cell cycle,” “Pathways in cancer,” and “p53 signaling pathway” (Figure 4B). In more details, eight risk lncRNAs (TYMSOS, ATP2A1-AS1, AC092718.4, MAFG-DT, AC108860.2, AC006329.1, LINC00665, AC099850.3, and AFAP1.AS1) from the “Cell cycle” pathway (Figure 4C) showed a ceRNA relationship with E2F1, 2, 3 (E2F1 and E2F2), Chk1 (CHEK1), Cdc25A (CDC25A), CycA (CCNA2), ORC (ORC1 and ORC2), and Dbf4 (DBF4) in SLGCeNBC, respectively. On the other hand, TYMSOS, AC092718.4, and AC006329.1 jointly upregulated the expression of CHEK1 (Chk1) and CCNA2 (CycA), inhibited the phosphorylation of Cdc25A (CDC25A), to reduce the dephosphorylation of CDK1 and CDK2, and further improved the phosphorylation levels of Rb and Dp-1,2. Consistently, ATP2A1-AS1, MAFG-DT, AC108860.2, and LINC00665 upregulated the expression of E2F1 and E2F2 (E2F1, 2, 3). To sum up, the expression of S-phase proteins CycE was indirectly co-promoted. On the other hand, the expression of DBF4 (Dbf4) was upregulated by TYMSOS, AC092718.4, and LINC00665, which accelerated the phosphorylation of MEM, and the ORC (origin recognition complex) was upregulated by AC092718.4, AC006329.1, AC099850.3, and AFAP1.AS1, which ultimately promoted DNA biosynthesis. For the “Pathway in cancer” pathway (Figure 4D), a total of 18 risk lncRNAs synergistically regulated ER (ESR1), E2F (E2F1 and E2F2), CyclinA1 (CCND1), and Survivin (BIRC5), which indirectly leads to tumor cell proliferation. In the eight risk lncRNAs, there are clear reports that LINC00665 is related to BRCA. LINC00665, which acted as ceRNA, promoted BRCA progression and induced an epithelial–mesenchymal transition-like phenotype via the competitively upregulation of LIN28B expression (Ding et al., 2020; Ji et al., 2020).

FIGURE 4

Figure 4. Synergistic regulation of the risk lncRNAs. (A) Risk lncRNA network; PCGs regulated by more than three risk lncRNAs were selected. The red and blue nodes represent lncRNAs and PCGs, respectively. (B) KEGG and GO enrichment analysis was performed using genes in the risk lncRNA network. The X-axis shows the –log10 of the P-value; the Y-axis shows the names of the pathway/term. The number at the top of the bar indicates the number of PCGs annotated to the pathway/term. (C) Risk lncRNAs regulate the cell cycle pathway (local). (D) Risk lncRNAs regulate the pathway in cancer (local).

By GO function analysis, we found that many biological processes were significant including “DNA replication,” “cell division,” “DNA repair,” “cell proliferation,” “DNA synthesis involved in DNA repair,” “regulation of cell cycle,” “DNA damage checkpoint,” “peptidyl-threonine phosphorylation,” “peptidyl-serine phosphorylation,” and “cytokinesis” (Figure 4B). Obviously, these biological processes are closely related to tumor cell proliferation differentiation and apoptosis. For the GO terms “peptidyl-threonine phosphorylation” and “peptidyl-serine phosphorylation,” Saeidi et al. (2020) have found that peptidyl-prolyl cis-trans isomerase NIMA-interacting 1 is highly overexpressed in human breast tumor tissues and H-Ras transformed human mammary epithelial (H-Ras MCF10A) and MDA-MB-231 BRCA cells.

Construction and Evaluation of Risk Prediction Model in the Training Dataset

We used the TCGA dataset to develop the risk prediction model and construct the prognostic signature, since we discovered 17 lncRNAs out of the 30 selected ones associated with the survival of BRCA patients. Then, we used the 17 prognostic lncRNAs to develop the risk prediction model and obtained 2¹⁷-1 = 131,071 risk models. We performed Kaplan–Meier analysis and compared the predictive ability of 131,071 signatures. A four-lncRNA signature (AP000851.1, LINC01977, MAFG-DT, and SIAH2-AS1) was found to have the minimum log-rank P-value (Figure 5A). The regression coefficients of the four lncRNAs (AP000851.1, LINC01977, MAFG-DT, and SIAH2-AS1) were all negative, which means they were related to BRCA poor prognosis (Supplementary Table 6).

FIGURE 5

Figure 5. Survival prediction performance of the four-lncRNA signature (AP000851.1, LINC01977, MAFG-DT, and SIAH2-AS1). (A) Kaplan–Meier curves for BRCA patients stratified by the four-lncRNA signature. P-values were calculated using a log-rank test. (B) ROC analysis to assess the survival accuracy of the four-lncRNA signature for 1, 3, 5, and 10 years of survival.

Survival Prediction Performance of the Four-lncRNA Signature in the TCGA Dataset

Each patient in the TCGA dataset received a risk score based on the four-lncRNA signature. Then, the patients with BRCA in the TCGA dataset were divided into high-risk (n = 537) or low-risk group (n = 538) based on the median risk score. Kaplan–Meier analysis demonstrated that patients in the low-risk group owned longer survival times than those in the high-risk group (median survival time: 11.69 vs. 9.48 years, log-rank test P < 0.001; Figure 5A). Subsequently, we performed time-dependent receiver operating characteristic (ROC) analysis to assess the survival accuracy of the four-lncRNA signature. In the TCGA dataset, the area under the ROC curve (AUC) for 1, 3, 5, and 10 years of survival were 0.68, 0.61 0.62, and 0.63, respectively (Figure 5B).

To test its prognostic independence, univariate and multivariable Cox regression analyses (Cox, 1972) were conducted. The multivariable Cox regression results in the TCGA datasets verified that the four-lncRNA signature can predict patients’ survival [high- vs. low-risk, hazard ratio (HR) training = 2.02, 95% confidence interval (CI) 1.43-2.86, P < 0.001, n = 1075; Table 3].

TABLE 3

Table 3. Cox regression analysis of the signature with BRCA survival (n = 1075).

Discussion

In recent years, the role of lncRNA has become a highly studied topic in the field of tumor research. Accumulating evidence indicates that lncRNA is involved in the oncogenesis and development of BRCA (Bray et al., 2018; DeSantis et al., 2019; Ji et al., 2020). With the development of high-throughput sequencing technology, a large number of lncRNA expression data involved in the occurrence and progression of cancer are emerging. It is not hard to infer that many of the lncRNAs related to BRCA may function as a complex and organized regulatory network for BRCA (Chi et al., 2019).

In this study, SLGCeNBC was constructed based on co-expression and miRNA-mediated RNA cross talks. Therefore, we identified 30 risk lncRNAs associated with BRCA by performing RWR in SLGCeNBC. In 30 risk lncRNAs, 13 risk lncRNAs have been confirmed to be associated with several cancers and nine risk lncRNAs have demonstrated a high association with BRCA. This shows that our method is effective and practical. By means of enrichment analysis using GO and KEGG, we found that these risk lncRNAs are significantly enriched in cancer-related biological processes and pathways, which have been stated by researchers (Hamilton and Infante, 2016; Anwar et al., 2018; Devanathan et al., 2020; Saeidi et al., 2020), and we found that the regulation of risk lncRNAs on BRCA-related genes is synergistic rather than alone. As a result, we found that these 30 risk lncRNAs are highly correlated with the occurrence and development of BRCA, and they may form the network system to jointly regulate the initiation and course of BRCA. These results suggest that risk lncRNAs may serve as novel diagnostic markers and treatment targets.

Importantly, we identified a four-lncRNA prognosis signature based on ceRNA network analysis, which could be used as the key clinical biomarker in BRCA prognosis. The poor prognosis of BRCA is mainly manifested in tumor metastasis, which is the leading cause of death (Zhang et al., 2017). Tumor metastasis is difficult to detect, which is discovered only when the tumor is large enough to be observed in regular follow-up imaging examination or to cause notable symptoms resulting from a tumor mass effect. This situation may exacerbate the patient prognosis. In our study, after survival analysis, we found that a four-lncRNA signature including four lncRNAs (AP000851.1, LINC01977, MAFG-DT, and SIAH2-AS1) has the ability to predict survival in patients with BRCA, which are expected to be novel predictors that may identify early metastasis. Additionally, we also investigate prognostic independence effects of the four-lncRNA signature in the TCGA dataset. As a result, the multivariable Cox regression results in the TCGA datasets verified that the four-lncRNA signature can predict patients’ survival. This result also demonstrated the clinical potential of the four-lncRNA signature.

Our methods also show some limitations. First of all, in this study, we conducted a bioinformatics analysis to identify the crucial factors in BRCA; results indicated that some genes (PCGs or lncRNAs) might play vital roles in the subtype cancers. Bioinformatics may infer only the functions of these lncRNAs; thus, it remains necessary to confirm the biological effects of these risk lncRNAs in BRCA in experimental studies. This result also encouraged us to validate the biological function and mechanism. In a further study, we will conduct the biological experiments to investigate these potential factors. Secondly, the risk lncRNAs identified here may not be all candidate risk lncRNAs associated with BRCA because of only limitations of lncRNA data. Finally, this study is based on the DE lncRNAs in BRCA. If all lncRNAs in tumors participate in the occurrence and development of tumors in the form of a network system, then some lncRNAs with no obvious changes may slip through the net and our results range may be narrowed. However, the research results of our study could contribute to accelerating the discovery of molecular biomarkers for diagnosis, treatment, and prognosis evaluation of BRCA.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: https://xenabrowser.net/datapages/.

Author Contributions

XS and LS contributed manuscript. LS and LW revised the manuscript. All authors contributed to the manuscript revision and, read and approved the submitted version.

Funding

This work was supported by the National Natural Science Foundation of China (31970636).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.664393/full#supplementary-material

Abbreviations

BRCA, breast cancer; lncRNA, long non-coding RNA; PCG, protein coding gene; ceRNA, competitive endogenous RNA; TCGA, The Cancer Genome Atlas; DE, differentially expressed; PCC, Pearson correlation coefficient; RWR, random walk with restart; SLGCeNBC, specific LncRNAs-PCGs ceRNA network for breast cancer; ROC, receiver operating characteristic; KM, Kaplan–Meier; AUC, area under the ROC curve; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; HR, hazard ratio; CI, confidence interval.

Footnotes

1. ^ https://genome-cancer.ucsc.edu
2. ^ www.r-project.org

References

Anwar, T., Arellano-Garcia, C., Ropa, J., Chen, Y. C., Kim, H. S., Yoon, E., et al. (2018). 38-mediated phosphorylation at T367 induces EZH2 cytoplasmic localization to promote breast cancer metastasis. Nat. Commun. 9:2801. doi: 10.1038/s41467-018-05078-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Bland, J. M., and Altman, D. G. (1998). Survival probabilities (the Kaplan-Meier method). BMJ 317:1572. doi: 10.1136/bmj.317.7172.1572

PubMed Abstract | CrossRef Full Text | Google Scholar

Bray, F., Ferlay, J., Soerjomataram, I., Siegel, R. L., Torre, L. A., and Jemal, A. (2018). Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J. Clin. 68, 394–424. doi: 10.3322/caac.21492

PubMed Abstract | CrossRef Full Text | Google Scholar

Chang, K. C., Diermeier, S. D., Yu, A. T., Brine, L. D., Russo, S., Bhatia, S., et al. (2020). MaTAR25 lncRNA regulates the tensin1 gene to impact breast cancer progression. Nat. Commun. 11:6438. doi: 10.1038/s41467-020-20207-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Chi, Y., Wang, D., Wang, J., Yu, W., and Yang, J. (2019). Long non-coding RNA in the pathogenesis of cancers. Cells 8:1015. doi: 10.3390/cells8091015

PubMed Abstract | CrossRef Full Text | Google Scholar

Chou, C. H., Shrestha, S., Yang, C. D., Chang, N. W., Lin, Y. L., Liao, K. W., et al. (2018). miRTarBase update 2018: a resource for experimentally validated microRNA-target interactions. Nucleic Acids Res. 46, D296–D302. doi: 10.1093/nar/gkx1067

PubMed Abstract | CrossRef Full Text | Google Scholar

Cox, D. R. (1972). Regression models and life tables. J. R. Statist. Assoc. 34, 187–220.

Google Scholar

Davis, J. E., Xie, X., Guo, J., Huang, W., Chu, W. M., Huang, S., et al. (2016). ARF1 promotes prostate tumorigenesis via targeting oncogenic MAPK signaling. Oncotarget 7, 39834–39845. doi: 10.18632/oncotarget.9405

PubMed Abstract | CrossRef Full Text | Google Scholar

DeSantis, C. E., Ma, J., Gaudet, M. M., Newman, L. A., Miller, K. D., Goding Sauer, A., et al. (2019). Breast cancer statistics, 2019. CA Cancer J. Clin. 69, 438–451. doi: 10.3322/caac.21583

PubMed Abstract | CrossRef Full Text | Google Scholar

Devanathan, N., Jones, S., Kaur, G., and Kimble-Hill, A. C. (2020). Using phosphatidylinositol phosphorylation as markers for hyperglycemic related breast cancer. Int. J. Mol. Sci. 21:2320. doi: 10.3390/ijms21072320

PubMed Abstract | CrossRef Full Text | Google Scholar

Di, J., Huang, H., Qu, D., Tang, J., Cao, W., Lu, Z., et al. (2015). Rap2B promotes proliferation, migration, and invasion of human breast cancer through calcium-related ERK1/2 signaling pathway. Sci. Rep. 5:12363. doi: 10.1038/srep12363

PubMed Abstract | CrossRef Full Text | Google Scholar

Ding, J., Zhao, J., Huan, L., Liu, Y., Qiao, Y., Wang, Z., et al. (2020). Inflammation-induced long intergenic noncoding RNA (LINC00665) increases malignancy through activating the double-stranded RNA-activated protein kinase/nuclear factor kappa b pathway in hepatocellular carcinoma. Hepatology 72, 1666–1681. doi: 10.1002/hep.31195

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, S., Hu, W., Zhao, Y., Zhou, H., Wen, W., Xu, M., et al. (2019). Long non-coding RNA MAGI2-AS3 inhibits breast cancer cell migration and invasion via sponging microRNA-374a. Cancer Biomark. 24, 269–277. doi: 10.3233/cbm-182216

PubMed Abstract | CrossRef Full Text | Google Scholar

Feng, C., Song, C., Ning, Z., Ai, B., Wang, Q., Xu, Y., et al. (2019). ce-Subpathway: identification of ceRNA-mediated subpathways via joint power of ceRNAs and pathway topologies. J. Cell Mol. Med. 23, 967–984. doi: 10.1111/jcmm.13997

PubMed Abstract | CrossRef Full Text | Google Scholar

Gibb, E. A., Brown, C. J., and Lam, W. L. (2011). The functional role of long non-coding RNA in human carcinomas. Mol. Cancer 10:38. doi: 10.1186/1476-4598-10-38

PubMed Abstract | CrossRef Full Text | Google Scholar

Goitre, L., Trapani, E., Trabalzini, L., and Retta, S. F. (2014). The Ras superfamily of small GTPases: the unlocked secrets. Methods Mol. Biol. 1120, 1–18. doi: 10.1007/978-1-62703-791-4_1

CrossRef Full Text | Google Scholar

Grimson, A., Farh, K. K., Johnston, W. K., Garrett-Engele, P., Lim, L. P., and Bartel, D. P. (2007). MicroRNA targeting specificity in mammals: determinants beyond seed pairing. Mol. Cell 27, 91–105. doi: 10.1016/j.molcel.2007.06.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, J. C., Li, C. Q., Wang, Q. Y., Zhao, J. M., Ding, J. Y., Li, E. M., et al. (2016). Protein-coding genes combined with long non-coding RNAs predict prognosis in esophageal squamous cell carcinoma patients as a novel clinical multi-dimensional signature. Mol. Biosyst. 12, 3467–3477. doi: 10.1039/c6mb00585c

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, J. C., Wu, Y., Chen, Y., Pan, F., Wu, Z. Y., Zhang, J. S., et al. (2018). Protein-coding genes combined with long noncoding RNA as a novel transcriptome molecular staging model to predict the survival of patients with esophageal squamous cell carcinoma. Cancer Commun. (Lond) 38:4. doi: 10.1186/s40880-018-0277-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Hamilton, E., and Infante, J. R. (2016). Targeting CDK4/6 in patients with cancer. Cancer Treat. Rev. 45, 129–138. doi: 10.1016/j.ctrv.2016.03.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Herrera-Solorio, A. M., Peralta-Arrieta, I., Armas López, L., Hernández-Cigala, N., Mendoza Milla, C., Ortiz Quintero, B., et al. (2020). LncRNA SOX2-OT regulates AKT/ERK and SOX2/GLI-1 expression, hinders therapy, and worsens clinical prognosis in malignant lung diseases. Mol. Oncol. 15, 1110–1129. doi: 10.1002/1878-0261.12875

PubMed Abstract | CrossRef Full Text | Google Scholar

Hothorn, T., and Zeileis, A. (2008). Generalized maximally selected statistics. Biometrics 64, 1263–1269. doi: 10.1111/j.1541-0420.2008.00995.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Jeggari, A., Marks, D. S., and Larsson, E. (2012). miRcode: a map of putative microRNA target sites in the long non-coding transcriptome. Bioinformatics 28, 2062–2063. doi: 10.1093/bioinformatics/bts344

PubMed Abstract | CrossRef Full Text | Google Scholar

Ji, W., Diao, Y. L., Qiu, Y. R., Ge, J., Cao, X. C., and Yu, Y. (2020). LINC00665 promotes breast cancer progression through regulation of the miR-379-5p/LIN28B axis. Cell Death Dis. 11:16. doi: 10.1038/s41419-019-2213-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Karagkouni, D., Paraskevopoulou, M. D., Tastsoglou, S., Skoufos, G., Karavangeli, A., Pierros, V., et al. (2020). DIANA-LncBase v3: indexing experimentally supported miRNA targets on non-coding transcripts. Nucleic Acids Res. 48, D101–D110. doi: 10.1093/nar/gkz1036

PubMed Abstract | CrossRef Full Text | Google Scholar

Labhart, P., Karmakar, S., Salicru, E. M., Egan, B. S., Alexiadis, V., O’Malley, B. W., et al. (2005). Identification of target genes in breast cancer cells directly regulated by the SRC-3/AIB1 coactivator. Proc. Natl. Acad. Sci. U.S.A. 102, 1339–1344. doi: 10.1073/pnas.0409578102

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J. H., Liu, S., Zhou, H., Qu, L. H., and Yang, J. H. (2014). starbase v2.0: decoding miRNA-ceRNA, miRNA-ncRNA and protein-RNA interaction networks from large-scale CLIP-Seq data. Nucleic Acids Res. 42, D92–D97. doi: 10.1093/nar/gkt1248

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, R., Peng, C., Zhang, X., Wu, Y., Pan, S., and Xiao, Y. (2017). Roles of Arf6 in cancer cell invasion, metastasis and proliferation. Life Sci. 182, 80–84. doi: 10.1016/j.lfs.2017.06.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Mao, Q., Lv, M., Li, L., Sun, Y., Liu, S., Shen, Y., et al. (2020). Long intergenic noncoding RNA 00641 inhibits breast cancer cell proliferation, migration, and invasion by sponging miR-194-5p. J. Cell Physiol. 235, 2668–2675. doi: 10.1002/jcp.29170

PubMed Abstract | CrossRef Full Text | Google Scholar

Meng, Q., Wang, X., Xue, T., Zhao, Q., Wang, W., and Zhao, K. (2020). Long noncoding RNA MIR99AHG promotes gastric cancer progression by inducing EMT and inhibiting apoptosis via miR577/FOXP1 axis. Cancer Cell Int. 20:414. doi: 10.1186/s12935-020-01510-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Nagini, S. (2017). Breast cancer: current molecular therapeutic targets and new players. Anticancer Agents Med. Chem. 17, 152–163. doi: 10.2174/1871520616666160502122724

PubMed Abstract | CrossRef Full Text | Google Scholar

Piñero, J., Ramírez-Anguita, J. M., Saüch-Pitarch, J., Ronzano, F., Centeno, E., Sanz, F., et al. (2020). The DisGeNET knowledge platform for disease genomics: 2019 update. Nucleic Acids Res. 48, D845–D855. doi: 10.1093/nar/gkz1021

PubMed Abstract | CrossRef Full Text | Google Scholar

Quinn, J. J., and Chang, H. Y. (2016). Unique features of long non-coding RNA biogenesis and function. Nat. Rev. Genet. 17, 47–62. doi: 10.1038/nrg.2015.10

PubMed Abstract | CrossRef Full Text | Google Scholar

Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 43:e47. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

Saeidi, S., Joo, S., Kim, S. J., Jagadeesh, A. S. V., and Surh, Y. J. (2020). Interaction between Peptidyl-prolyl Cis-trans isomerase NIMA-interacting 1 and GTP-H-Ras: implications for aggressiveness of human mammary epithelial cells and drug resistance. J. Cancer Prev. 25, 234–243. doi: 10.15430/jcp.2020.25.4.234

PubMed Abstract | CrossRef Full Text | Google Scholar

Soudyab, M., Iranpour, M., and Ghafouri-Fard, S. (2016). The role of long non-coding RNAs in breast cancer. Arch. Iran Med. 19, 508–517.

Google Scholar

Trovero, M. F., Rodríguez-Casuriaga, R., Romeo, C., Santiñaque, F. F., François, M., Folle, G. A., et al. (2020). Revealing stage-specific expression patterns of long noncoding RNAs along mouse spermatogenesis. RNA Biol. 17, 350–365. doi: 10.1080/15476286.2019.1700332

PubMed Abstract | CrossRef Full Text | Google Scholar

Ulitsky, I., and Bartel, D. P. (2013). lincRNAs: genomics, evolution, and mechanisms. Cell 154, 26–46. doi: 10.1016/j.cell.2013.06.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, W., Lou, W., Ding, B., Yang, B., Lu, H., Kong, Q., et al. (2019). A novel mRNA-miRNA-lncRNA competing endogenous RNA triple sub-network associated with prognosis of pancreatic cancer. Aging (Albany NY) 11, 2610–2627. doi: 10.18632/aging.101933

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiao, M., Liu, J., Xiang, L., Zhao, K., He, D., Zeng, Q., et al. (2020). MAFG-AS1 promotes tumor progression via regulation of the HuR/TB1 axis in bladder urothelial carcinoma. Clin. Transl. Med. 10:e241. doi: 10.1002/ctm2.241

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, X., Tang, S. C., Cai, Y., Pi, W., Deng, L., Wu, G., et al. (2016). Suppression of breast cancer metastasis through the inactivation of ADP-ribosylation factor 1. Oncotarget 7, 58111–58120. doi: 10.18632/oncotarget.11185

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, X., Yuan, X., Ni, J., Guo, J., Gao, Y., Yin, W., et al. (2021). MAGI2-AS3 inhibits breast cancer by downregulating DNA methylation of MAGI2. J. Cell Physiol. 236, 1116–1130. doi: 10.1002/jcp.29922

PubMed Abstract | CrossRef Full Text | Google Scholar

Yin, C., Mou, Q., Pan, X., Zhang, G., Li, H., and Sun, Y. (2018). MiR-577 suppresses epithelial-mesenchymal transition and metastasis of breast cancer by targeting Rab25. Thorac. Cancer 9, 472–479. doi: 10.1111/1759-7714.12612

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, J., Feng, C., Song, C., Ai, B., Bai, X., Liu, Y., et al. (2018). Identification and analysis of a key long non-coding RNAs (lncRNAs)-associated module reveal functional lncRNAs in cardiac hypertrophy. J. Cell Mol. Med. 22, 892–903. doi: 10.1111/jcmm.13376

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Wang, D. L., Yan, H. Y., Liao, J. Y., He, J. H., Hu, K. S., et al. (2017). Genome-wide study of ER-regulated lncRNAs shows AP000439.3 may function as a key regulator of cell cycle in breast cancer. Oncol. Rep. 38, 3227–3237. doi: 10.3892/or.2017.5975

PubMed Abstract | CrossRef Full Text | Google Scholar