|Home | About | Journals | Submit | Contact Us | Français|
Triple negative breast cancer (TNBC) is a particular subtype of breast malignant tumor with poorer prognosis than other molecular subtypes. Currently, there is increasing focus on long non-coding RNAs (lncRNAs), which can act as competing endogenous RNAs (ceR-NAs) and suppress miRNA functions involved in post-transcriptional regulatory networks in the tumor. Therefore, to investigate specific mechanisms of TNBC carcinogenesis and improve treatment efficiency, we comprehensively integrated expression profiles, including data on mRNAs, lncRNAs and miRNAs obtained from 116 TNBC tissues and 11 normal tissues from The Cancer Genome Atlas. As a result, we selected the threshold with |log2FC|>2.0 and an adjusted p-value >0.05 to obtain the differentially expressed mRNAs, miRNAs and lncRNAs. Hereafter, weighted gene co-expression network analysis was performed to identify the expression characteristics of dysregulated genes. We obtained five co-expression modules and related clinical feature. By means of correlating gene modules with protein–protein interaction network analysis that had identified 22 hub mRNAs which could as hub target genes. Eleven key dysregulated differentially expressed micro RNAs (DEmiRNAs) were identified that were significantly associated with the 22 hub potential target genes. Moreover, we found that 14 key differentially expressed lncRNAs could interact with the key DEmiRNAs. Then, the ceRNA crosstalk network of TNBC was constructed by utilizing key lncRNAs, key miRNAs, and hub mRNAs in Cytoscape software. We analyzed and described the potential characteristics of biological function and pathological roles of the TNBC ceRNA co-regulatory network; also, the survival analysis was performed for each molecule. These findings revealed that ceRNA crosstalk network could play an important role in the development and progression for TNBC. In addition, we also identified that some molecules in the ceRNA network possess clinical correlation and prognosis.
Breast cancer is a heterogeneous disease with varying histological abnormalities, cytogenetic aberrations, responses to therapy, and prognoses, and accounts for 16% of all female cancers and 22.9% of invasive cancers in women.1 Triple negative breast cancer (TNBC), a particular subtype of BC, is tested with negative for estrogen receptors (ER−), progesterone receptors (PR−), and Erb-B2 receptor tyrosine kinase 2 (HER2−).2 This kind of breast cancer generally has overall poor outcomes, such as more aggressive, higher distant recurrence and poorer survival in metastatic disease.3 Moreover, current targeted therapy strategy is available for most breast cancers that are ER or PR positive, but there are no specific targeted therapy options available for TNBC.4 Therefore, it is urgent and necessary that new treatments are found based on the specific mechanisms of TNBC to improve treatment efficiency and avoid the adverse effects of conventional methods.
Advances in next-generation RNA sequencing technology have brought to light the complexity of our genome. Non-coding RNAs (ncRNAs), such as micro RNAs (miRNAs) and the long non-coding RNAs (lncRNAs), accounting for the majority (98%) of the transcriptome, are defined as gene transcripts with important biological functions.5 Investigating the world of RNA is one of the most important challenges facing biologists today, and lncRNAs can represent potential new biomarkers and drug targets. Salmena et al’s competing endogenous RNA (ceRNA) hypothesis about the multiple molecular involved in post-transcriptional regulatory networks has received much attention in recent years; their study indicated that lncRNAs, cirRNAs, and pseudogenes can act as ceRNAs to suppress miRNA functions and communicate with mRNAs by harboring one or more miRNA response elements (MREs).6 An increasing number of experimental reports are supporting this advanced hypothesis. One such example is lnc-MD1, a muscle-specific lncRNA that has been shown to control muscle differentiation by acting as miR-135 and miR-133 “sponge,” and thus regulating the expression of MAML1 and MEF2C.7 In addition, Wang et al has identified Linc-ROR, which could serve as a ceRNA and acts as miR-145 “sponge”, thereby regulating SOX2, OCT4 and NANOG expression to control self-renewal and maintain pluripotency of human embryonic stem cells.8 However, to date, there is still little research on ceRNA-related lncRNAs’ potential roles in relevant regulatory mechanisms of TNBC.
The Cancer Genome Atlas (TCGA) project was launched in 2005 and has accelerated the comprehensive understanding of cancer genomic profiles that has improved diagnostic methods, therapy standards, and preventive strategies. By the end of 2015, TCGA had expanded to multiple types, reaching 30 different tumor types with genome expression profiles and clinical data through the large-scale public “atlas”.9 Also, the latest progress in the integrated multi-dimensional analyses development of technology bioinformatics has shed new insight on the cancer profiling. Weighted gene co-expression network analysis (WGCNA) is a systems biology algorithm for constructing the correlation patterns among genes across samples.10 Furthermore, WGCNA can be used for uncovering and clustering highly correlated molecules into separate modules that provide information on hub nodes and external sample traits with related modules. The varied studies of WGCNA are widely used in cancer research, genetics of species and other complex diseases to investigate and identify related genes, biological processes, and candidate biomarkers. However, few studies have utilized the open access TCGA database for cancer research with WGCNA method. This research has mined TNBC data in the TCGA database that analyzed the integrated RNA expression profiles from 116 TNBC samples and 11 normal samples. We used the WGCNA and protein–protein interaction (PPI) network to identify hub mRNAs from the differentially expressed mRNAs (DEmRNAs) profile. The databases, miRTarBase and miRcode, were used to screen the key differentially expressed lncRNAs (DElncRNAs) and differentially expressed micro RNAs (DEmiRNAs), as well as DEmiRNAs-DEmRNAs interactions. Hereafter, we constructed a ceRNA network to elaborate on the interactions and potential crosstalk between the hub DElncRNAs, DEmiRNAs and DEmRNAs. In addition, relevant survival analyses were performed to determine the prognostic genes that possess clinical traits.
All RNA sequencing data sets were included in this study that met the following criteria: 1) detailed clinical data consisted of age, gender, race, TNM stage, and histological subtype and grade; 2) patients possessed a certain follow-up time. Moreover, details of set grouping conditions were as follows: the tumor group included the breast cancer samples of triple negative; and the control group included normal physical tissue. Overall, publicly available breast cancer level 3 RNAseq and miRNASeq data derived from 127 samples, including 116 TNBC tumor tissues and 11 adjacent normal tissues, were obtained from TCGA. These data were classified into two cohorts. According to the TCGA ethics committee, the data and analysis of this present study did not require further approval.
We merged tumor sample and normal sample data and deleted expressed data that closed to zero. The DEmRNAs and DEmiRNAs of raw count data were processed through the limma package, which assessed the data using linear models and empirical Bayesian methods.11 The threshold of differential expression was defined with an absolute log 2-fold change (FC) >2 and an adjusted p-value <0.05.
In addition, the GENCODE lncRNA annotation (V22) was used to discover and define the lncRNAs from DERNAs.12 In our study, lncRNAs not included in the annotation were discarded.
The differential expression genes were applied to find the scale-free gene modules of co-expression and highly correlated genes constructed by WGCNA package.10 The outlier samples were identified and removed to ensure the reliability of co-expression network construction. To build the scale-free co-expression relationship between gene pairs, an appropriate soft threshold power parameter was set to 4, which could penalize weaker connections and strengthen stronger connections. Topological overlap matrix (TOM) was performed by adjacency transformation. Then, we calculated 1-TOM as the distance measured and through the dynamic tree cut method to identify hierarchical clustering genes and modules. The minimum size cut-off of 100 for the resulting dendrogram of each module, and below which separate modules would be merged, was set to 0.25.
In addition, for investigating the correlation between gene expression modules and clinical features, we plotted the heatmap of all modules with the relation cut-off and p-value at eight different conditions, which included pathology diagnosis age, laterality, histological type, positive lymph nodes, pathology stage, pathology T stage, pathology M stage, pathology N stage. The previously mentioned data were analyzed by Pearson’s correlation tests.
In order to further explore the characteristics of each module, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of differentially expressed genes (DEGs) were identified by KEGG Orthology Based Annotation System (KOBAS) web server, which provides abundant annotation information of gene sets from multiple databases covering pathways (KEGG PATHWAY, Reactome, Biocyc, Panther) and supports gene and protein functional annotation.13,14 In addition, the significant enrichment KEGG pathway should have contained at least three DEGs and had a false discovery rate (FDR) of <0.05. We all know that protein interaction data analysis serves as an important tool to unravel the molecular basis of a disease. Thus, it is feasible to discover TNBC-related DEGs of different modules based on protein network. This study constructed interaction networks and retrieved data from STRING (Search Tool for the Retrieval of Interacting Genes/Proteins, http://string.embl.de/), a well-known online database integrating known and predicted protein interactions.15 Additionally, we set the threshold for interaction score at >0.9, which indicates the highest confidence of two protein interaction occurrences.
The hub dysregulated genes and ncRNA associated ceRNA crosstalk network in TNBC was constructed as follows: 1) prediction of target mRNAs of hub mRNAs: we selected the hub DEmRNAs from WGCNA modules as the target genes, and found out the miRNA-mRNA pairs with experimental support by miRTarBase (http://mirtarbase.mbc.nctu.edu.tw/).16 Then, hub DEmRNAs were mapped into the interactions. 2) Prediction of target lncRNAs of miRNAs: the lncRNA-miRNA interactions were predicted through miRcode (http://www.mircode.org/) database.17 The DElncRNAs and DEmiRNAs were mapped into the interactions. 3) Construction of ceRNA network: TNBC-specific lncRNA-miRNA-mRNA interactions were revealed. The ceRNA crosstalk network was constructed with the previously mentioned interactions. The construction of the network was visualized using the Cytoscape 3.4.0 software.
Furthermore, to explore potential biological mechanisms of the predicted ceRNA crosstalk network, the KOBAS web server was used to annotate KEGG pathways of hub genes that could be activated by key lncRNAs and miRNAs. Each item of enrichment had a threshold of FDR <0.05. The network of pathways was visualized using the Cytoscape software.
To identify prognostic lncRNA and miRNA signatures, combining the clinical data of those patients in TCGA, we established the life curves of those samples with each node in the ceRNA network. Survival analysis was performed using a standard Kaplan–Meier univariate curve through the “survival” package in R 184.108.40.206 Molecule with p-value <0.05 was considered as having statistical significance.
The level three expression profiles of TNBC and corresponding clinical data (retrieved up to July 2017) were downloaded by Data Transfer Tool (https://tcga-data.nci.nih.gov/tcga/) in the TCGA database. Data of 127 patients were available. Moreover, to obtain DEmRNAs and DEmiRNAs, we applied the limma algorithm, which can identify the different expressions of TNBC and normal samples. Finally, a total of 2,791 DEmRNAs were identified by the threshold of log2FC >2 and p-value <0.01, of which 1,419 were down-regulated and 1,372 were up-regulated. In the same way, 67 DEmiRNAs were identified by limma, of which 46 were up-regulated in cancer samples and 21 down-regulated.
The GENCODE lncRNA annotation (V22) was used to discover and define the lncRNAs from DERNAs. According to GENECARDS annotation, we identified a total of 723 aberrantly expressed lncRNAs in TNBC tissues compared to the normal tissues. Two hundred and eighty-eight (39.83%) of the DElncRNAs were up-regulated, and 435 (60.17%) were down-regulated. We have shown the top 10 mRNAs, miRNAs and lncRNAs with |log2FC|>2.0 in Table 1.
To find the functional clusters in TNBC samples, we took the different regulated mRNAs data to build the correlation network through WGCNA. The hierarchical dendrogram of WGCNA modules, with different colors representing different modules, was constructed with an appropriate soft threshold (Figure 1). Eventually, it gained a total of five color modules, including red, blue, yellow, green, and brown. Additionally, the region of color gray represented those genes that were not divided into modules.
To further understand the characteristics of each module, the annotation of KEGG pathways for each module is listed in Table S1. Furthermore, we constructed the PPI networks to screen out hub DEmRNAs. As is shown in Figure 2, there are 25 nodes and 35 interactions exhibited in the red module. According to the highest confidence of 0.9 as threshold, two hub genes, EFNA3, and FOXP3, were noted in the center in the PPI network. In addition, we found that EFNA3 was mediated by four pathways, including “MicroRNAs in cancer,” “PI3K-Akt signaling pathway,” “Rap1 signaling pathway,” and “Ras signaling pathway.” In the blue module, we found 99 nodes and 869 interactions. And, four hub genes were significantly enriched in the blue PPI network, including CDC25A, CCNE2, STMN1 and CCNA2. The KEGG pathway annotation of blue module revealed the enrichment of genes involved in “cell cycle” (CDC25A, CCNE2, CCNA2), “MicroRNAs in cancer” (CDC25A, CCNE2) and “p53 signaling pathway” (CCNE2).
The yellow module contained 21 nodes and 56 interactions. Also, three hub genes which are members of the peptidase M10 family of matrix metalloproteinases (MMPs), including MMP1, MMP9 and MMP13, were noticed in the PPI network. We identified the KEGG pathways of this module related to critical cellular processes and organismal systems (endocrine system). In the green module, we obtained 125 nodes and 451 interactions. A total of 10 hub genes with significant centrality values were identified in this module, consisting of ABRACL, ABCB11, CHEK1, CD34, FOS, IL6, IFNB1, ESR1, FGF2 and SAMD5. Pathway enrichment analysis of the green module revealed that the enriched terms focused on environment signaling progress and the endocrine system, like “Neuroactive ligand-receptor interaction,” “cAMP signaling pathway,” and “Regulation of lipolysis in adipocytes.” There were 39 nodes and 32 interactions in the brown module. We identified four nodes in center that serve as hub genes in the PPI network, which are as follows: PTHLH, HMGA2 and AKAP12. Moreover, in this module, pathway enriched terms were concentrated in metabolism progress and disease-related pathways; for instance, the hub gene HMGA2 was enriched in “MicroRNAs in cancer”, “Transcriptional misregulation in cancer”.
After gaining the summarized clinical and pathological information (Table 2), we constructed the feature heatmap for investigating the relation between modules and clinical characteristics (Figure 3). As a result, we noticed that the red module for two variables, encompassing positive lymph nodes (r=0.49, p=0.002) and pathology M stage (r=0.54, p=0.009), displayed the highest correlation coefficient. Moreover, the yellow module also showed the higher association in pathology M stage (r=0.34), though no statistical significance was revealed (p=0.4). The green module exhibited the highest correlation in pathology N stage (r=0.24, p=0.01). Taken together, the DEmRNAs of the previously mentioned modules may be the key regulators for the development of TNBC.
To further investigate the role and relation of lncRNA and miRNA mediation in TNBC, we searched for DEmiRNAs targeted by these 22 hub mRNAs obtained from WGCNA and PPI network analysis. Through the predicted miRNAs and algorithm analyzed dysregulated expressed 67 miRNAs. As a result, 11 intersection miRNAs related to the 22 hub genes were identified by miRTarBase database, and some of them have been reported to be cancer-associated miRNAs such as hsa-mir-145, hsa-mir-9, and hsa-mir-503. The miRNAs target with lncRNA through MREs in the ceRNA network, and thus, we detected whether key miRNAs and the above 723 differently expressed lncRNAs have the potential MREs based on the information retrieved from the miRcode database. The predicted results showed that the key miRNAs interacted with 14 DElncRNAs, which could act as key lncRNAs. Finally, we obtained 22 hub DEmRNAs mediating with 11 DEmiRNAs and 14 DElncRNAs, which could be the key nodes involved in the ceRNA network.
Based on combination of the above data, the lncRNA-miRNA-mRNA ceRNA network was established and visualized using the Cytoscape software (Figure 4A). The detailed information of each node is shown in Table 3. In the ceRNA network, some of them have been reported to be cancer-associated moleculars such as MMP1, ESR1, STMN1 and HOTAIR. We noticed that DGCR5 is the key lncRNA interrelated with seven key miRNAs (hsa-mir-145, hsa-mir-33a, hsa-mir-142, hsa-mir-144, hsa-mir-503, hsa-mir-183, and hsa-mir-210), and may serve as ceRNA to inhibit target miRNAs, and mediated related hub genes translation such as HMGA2, CCNA2, SAMD5, ABRACL, and IFNB1. And, HNF1A-AS1 may act as an endogenous sink by binding hsa-mir-183, hsa-mir-141, hsa-mir-145 and may be associated with targeted genes transcripts such as AKAP12, ESR1, MMP1, CCNA2. In addition, RP11-576D8.4 is a lncRNA, with no current research to find its function. In our study, it may have played a role in TNBC by its interactions with hsa-mir-145, hsa-mir-210, which then regulate the target hub mRNAs FOXP3, STMN1, EFNA3, ESR1, MMP1, and IFNB1. Moreover, EMX2OS, LINC00461, and LINC00484 also serve as the center molecule of several miRNAs. Thus, we suppose that DGCR5, HNF1A-AS1, LINC00484, LINC00461, RP11-576D8.4 and EMX2OS may closely contribute to the pathogenesis of TNBC.
To evaluate whether the ceRNA network involved in specific enrichment biological pathways could be activated by the 19 key DEmiRNAs and DElncRNAs, we annotated by KOBAS web server, a total of significantly enriched 14 KEGG pathways (FDR <0.05) (Table 4). Then, a network involved in key DElncRNAs, DEmiRNAs, and signal pathways interactions was constructed using Cytoscape (Figure 4B). Several critical tumor-related signaling pathways were involved in our ceRNA crosstalk network; for instance, “Pathways in cancer, MicroRNAs in cancer, p53 signaling pathway and PI3K-Akt signaling pathway”. Also, revealed were the enriched terms focused on critical cellular processes and organismal systems, such as “cell cycle, Estrogen signaling pathway, endocrine resistance, EGFR tyrosine kinase inhibitor resistance”.
Among these key DElncRNAs, we noticed that three lncRNAs had more interactions than other regulators, and they may be key regulators related to the progression of TNBC, including DGCR5, RP11-576D8.4, and EMX2OS. It is interesting that DGCR5 has enriched a total of 21 pathogenesis-related signaling pathways by mediating six DEmiRNAs, such as “MAPK signaling pathway”, “EGFR tyrosine kinase inhibitor resistance”, and “Transcriptional misregulation in cancer”. The lncRNA EMX2OS can directly regulate hsa-mir-210 and hsa-mir-503 to activate 12 pivotal pathways, such as “p53 signaling pathway”, “Proteoglycans in cancer”, and “Ras signaling pathway”. In addition, we found that RP11-576D8.4, a novel lncRNA, can mediate hsa-mir-145 and hsa-mir-210 to regulate 19 pathways, including “Endocrine resistance”, “Estrogen signaling pathway”, “MicroRNAs in cancer”.
The Kaplan–Meier curves were performed to correlate expression pattern and clinical information of all TNBC samples to identify key nodes significantly associated with patients’ survival (p<0.05). A total of 22 hub genes, 14 key lncRNAs and 12 key miRNAs involved in ceRNA network were subjected to survival analysis. Finally, we noticed that the over-expression levels of seven key molecules (AKAP12, FOS, EMX2OS, MYCNOS, RP11-542B15.1, hsa-mir-9, and hsa-mir-183) predicted patients with poor outcome. On the contrary, the patients with high expression of hsa-mir-145, LINC00461, RP11-576D8.4 and RP11-496D24.2 have shown longer overall survival time than patients with the lower expression. All survival analyses plotted are shown in Figures 5 and and66.
TNBC is a subtype of breast malignant tumor with poorer prognosis than other molecular subtypes. Even the development of chemotherapy has made progress and the patients can achieve pathological complete response after neoadjuvant chemotherapy, TNBC patients still have a higher metastasis rate and poorer prognosis than other breast cancer subtypes.19 With the increasing amounts of evidence indicating that ncRNAs have important roles in various biological processes, the correlated ceRNA networks have also been implicated in many kinds of tumor and regarded as a mechanism potentially driving tumors progression and metastasis,20,21 including in lung cancer,22 gastric cancer,23 and hepatocellular carcinoma.24 In addition, the ceRNA crosstalk could also relate to the exact regulatory mechanism of TNBC. Therefore, to find the proven effective target therapies for TNBC, we comprehensively integrated mRNA and miRNA data from TCGA, adopting WGCNA and other methods to identify the lncRNA-miRNA-mRNA interactions mediated in the ceRNA network for further investigation of the regulatory mechanisms of lncRNAs in this study.
lncRNAs are defined as non-protein coding transcripts longer than 200 nucleotides. Nevertheless, much research in recent years has shown that some lncRNAs can encode proteins to regulate cellular processes. In our work, 11 hub lncRNAs were identified as hub lncRNAs, and six of them (DGCR5, LINC00484, LINC00461, EMX2OS, RP11-576D8.4, and HNF1A-AS1) were distinctly enriched in the construction of the ceRNA network by integrated analysis. After detailed observation of the ceRNA network, we noticed that DGCR5 is the center of the network that may be involved in regulating the target genes by competing for five common shared miRNAs (hsa-mir-145, hsa-mir-141, hsa-mir-144, hsa-mir-183, and hsa-mir-210). In previous studies, DGCR5 was shown to deregulate expression in various tumors that could serve as a potential biomarker for diagnosis and prognosis. Sun et al had identified that DGCR5 was also down-regulated in pancreatic ductal adenocarcinoma cell lines, and its dysregulated expression could interact with miR-320a to inhibit tumor cell proliferation and migration.25 Huang et al had found that down-regulation of lncRNA DGCR5 correlated with poor prognosis in hepatocellular carcinoma.26 In our study, we noticed that DGCR5 with high-expression can compete with miRNAs to regulate the expression of the WGCNA green and red module target genes, such as ESR1, CHEK1, SAMD5, FGF2, EFNA3 and FOXP3 involved in the ceRNA network. Moreover, combined with clinical characteristics, the WGCNA green and red module demonstrated that they had the highest correlation coefficient in pathology N stage, pathology M stage, and positive lymph nodes of TNBC. Thus, it was hypothesized that knocking down the expression of DGCR5 may inhibit cell proliferation and lymphatic metastasis in TNBC. However, the survival analysis for DGCR5 showed no significant p-value, and these findings require further research to annotate its impact on breast cancer.
Besides, LINC00461 has the second-largest interactions in the ceRNA network. Recent study identified that dysfunctional expression of LINC00461 was significantly upregulated in both human and mouse glioma tissues. Deguchi et al demonstrated that LINC00461 could be a potential oncogene to regulate topoisomerase 2 alpha with sponging miR-411-5p in glioma.27 In this study, we identified that LINC00461 with high-expression may compete with three hub DEmiRNAs (miR-503, miR-144, miR-141) to regulate the expression of the target genes. In addition, the survival analysis showed that patients with highly expressed LINC00461 have prolonged overall survival time. In the construction of the ceRNA network, another two key DElncRNAs, RP11-496D24.2 and RP11576D8.4, were also related to the clinical outcome of TNBC patients. Low expression of these two lncRNAs suggested a poor prognosis in patients. It can be hypothesized that inducing the expression of the lncRNAs, including LINC00461, RP11-496D24.2 and RP11576D8.4, can prevent cell proliferation, and they could become the potential diagnosis and treatment targets for TNBC.
Even though miRNAs are defined as small ncRNA molecules, which have been well investigated in various human diseases, only a few lncRNAs have been functionally and mechanically characterized. In this study, we constructed the new ceRNA network to comprehensively illuminate the relations of lncRNAs and miRNAs. In our study, we identified 11 miRNAs, namely, mir-9, mir-33a, mir-139, mir-141, mir-142, mir-144, mir-145, mir-183, mir-210, mir-503 and let-7c, which should be further investigated. The mir-9, mir-183 and mir-145 were particularly important in the core of the ceRNA network. Tang et al verified that exosome mir-9 can suppress the expression of HBGF-5 and be a functional biomarker in hepatocellular carcinoma.28 Our study revealed that over-expression of mir-9 predicted a shorter lifetime compared to patients with lower mir-9 expression in TNBC. As for mir-183, Peng et al showed that mir-183 can be a potential biomarker for lung adenocarcinoma.29 Our survival analysis indicated that high expression of mir-183 was negatively linked to overall survival of TNBC patients. Regarding mir-145, Wang et al proved that the loss of mir-145 can induce Ras signaling pathway to promote aggressive Pten-deficient basal-like breast cancer.30 Our data also confirmed that the low expression of mir-145 predicted a shorter lifetime in TNBC. However, these findings need more research to investigate whether these aberrant miRNAs have a specific role in TNBC tumorigenesis. Many protein coding genes involved in the new ceRNA network were recently reported as key oncogenes participating in breast cancer development and progression, such as CCNE2, ESR1, HMGA2 and STMN1, which have also proven promising therapeutic targets in TNBC.31–33 Moreover, another two DEmRNAs, namely, FOS and AKAP12, demonstrated interactions with hub DEmiRNAs, and their high-expression predicted poor outcome in TNBC.
Based on the results of co-expression by KEGG pathways network, we could identify how the lncRNAs from ceRNA network activated the crucial pathways to regulate TNBC development and progression. Two enriched pathways, PI3K-Akt signaling pathway and mitogen-activated protein kinase (MAPK) signaling pathway, were reported as critical pathways in TNBC. Previous studies have identified that the crosstalk between the PI3K/Akt and MAPK pathways can overly activate to cause cell proliferation/metabolism, which impacts the prognosis of TNBC patients.34 Among the KEGG pathways network, five DElncRNAs, including DGCR5, EMX2OS, LINC00461, LINC00484, and RP11-576D8.4, may play specific roles for tumor-suppressing effects in TNBC, via down-regulation of two signal pathways.
Currently, although lncRNA-miRNA-mRNA network is receiving increasingly more attention, investigations focusing on that the way remarkable ncRNA mediate human diseases. Chen et al constructed a lncRNA-miRNA-mRNA ceRNA network to provide new insight into the common and specific subtype ceRNA interactions, which investigate the different categories of breast tumor (luminal A, luminal B, HER-2 enriched and TNBC subtype).35 However, this study only focused on different ncRNA data, and there were no combinations with clinical data to predict the prognosis in TNBC. Additionally, a similar study reported by Zhang et al also used the other kind of network-based method to analyse four different TCGA datasets for TNBC.36 The research used different molecular features, including expression profiles, gene copy number variants datasets and DNA methylation datasets, to identify some molecules that may serve as potential biomarkers.37 Nonetheless, the previous reports still could not provide an integrated WGCNA method to construct the ceRNA network for TNBC. And thus, the lncRNAs, miR-NAs, or other types of biological molecules have not been mechanically and functionally characterized.
Our research has some limitations. After a comprehensive search of Gene Expression Omnibus database, an international public repository that archives and freely distributes microarray, next-generation sequencing, and other forms of high-throughput functional genomic datasets, we found that there is no such independent data set that contains both coding and non-coding expression profiles for TNBC. Therefore, we lacked the validation step using independent datasets to demonstrate our findings. Moreover, since our research was only based on limited samples and information from TCGA databases, there is no more detailed information in TCGA for considering the subtypes analysis. Thus, further research is needed for verifying our study and unraveling the mechanism of ceRNA in TNBC.
Taken together, our study used an integrative biology approach to analyze non-coding and coding RNA expression profiles between normal and tumor samples from the TCGA database. Moreover, we performed WGCNA algorithm to select hub genes and construct the new lncRNA-miRNA-mRNA ceRNA crosstalk network to display a novel regulatory mechanism for further investigation of TNBC.
We gratefully acknowledge The Cancer Genome Atlas pilot project for providing all data in this research. Additionally, we are very grateful to Dr Chuiguo Huang for constructive and insightful suggestions. Financially supported by National Natural Science Foundation of China (nos 81673979, 81473688, 81373314, 81173265); Education Program of China for New Century Excellent Talents (no NCET-13-0827); Traditional Chinese Medicine Administration Project of Guangdong Province, China (no 20141070); Science and Technology Program of Guangzhou, China (no 2014J4100104); Science and Technology Program of Guangdong, China (nos 2014A020212672, and 2014A020210015); Natural Science Foundation of Guangdong Province, China (nos 2016A030313114; 2015A030313333); National Undergraduate Training Programs for Innovation and Entrepreneurship in 2015 (no 201510559046); Guangzhou Municipal Enterprise Research and Development Institutions Construction Project (no 201503010064); The Fundamental Research Funds for the Central Universities\Research and Cultivation and Innovation Fund of Jinan University, Guangzhou, China (nos 21615412, 21615464, and 21617467); Guangzhou Municipal Science and Technology Program in 2018 (project leader: Min Ma).
The authors report no conflicts of interest in this work.