Identification of potential key genes and pathways predicting pathogenesis and prognosis for triple-negative breast cancer

Background Triple negative breast cancer (TNBC) is a specific subtype of breast cancer with a poor prognosis due to its aggressive biological behaviour and lack of therapeutic targets. We aimed to explore some novel genes and pathways related to TNBC prognosis through bioinformatics methods as well as potential initiation and progression mechanisms. Methods Breast cancer mRNA data were obtained from The Cancer Genome Atlas database (TCGA). Differential expression analysis of cancer and adjacent cancer, as well as, triple negative breast cancer and non-triple negative breast cancer were performed using R software. The key genes related to the pathogenesis were identified by functional and pathway enrichment analysis and protein–protein interaction network analysis. Based on univariate and multivariate Cox proportional hazards model analyses, a gene signature was established to predict overall survival. Receiver operating characteristic curve was used to evaluate the prognostic performance of our model. Results Based on mRNA expression profiling of breast cancer patients from the TCGA database, 755 differentially expressed overlapping mRNAs were detected between TNBC/non-TNBC samples and normal tissue. We found eight hub genes associated with the cell cycle pathway highly expressed in TNBC. Additionally, a novel six-gene (TMEM252, PRB2, SMCO1, IVL, SMR3B and COL9A3) signature from the 755 differentially expressed mRNAs was constructed and significantly associated with prognosis as an independent prognostic signature. TNBC patients with high-risk scores based on the expression of the 6-mRNAs had significantly shorter survival times compared to patients with low-risk scores (P < 0.0001). Conclusions The eight hub genes we identified might be tightly correlated with TNBC pathogenesis. The 6-mRNA signature established might act as an independent biomarker with a potentially good performance in predicting overall survival. Electronic supplementary material The online version of this article (10.1186/s12935-019-0884-0) contains supplementary material, which is available to authorized users.


Background
Triple-negative breast cancer (TNBC) is defined as a subtype of aggressive breast cancer, accounting for 10-20% of all breast cancer cases [1]. TNBC subjects lack expression of the estrogen receptor (ER) and progesterone receptor (PR) and does not amplify the human epidermal growth factor receptor 2 (HER2) [2]. TNBC is more commonly diagnosed among young women and is more prone to relapse and visceral metastasis, compared with other breast cancer subtypes [3][4][5]. Due to the absence of molecular targets, patients diagnosed with TNBC cannot receive endocrine or HER2 targeted therapy [6], increasing the difficulty of treatment for them [7]. Chemotherapy is still the main adjuvant treatment option for patients with TNBC [8]. TNBC remains a disease associated with poor prognosis and limited treatment options because many tumours are resistant to chemotherapy and rapidly relapse or metastasize after adjuvant therapy [9]. The identification of uniform targets can help achieve more effective and less toxic treatment. Hence, it is imperative and urgent to explore new therapeutic targets for TNBC [10].
Recently, many biomarkers have been developed for breast cancer. For example, CD82, a potential diagnostic biomarker for breast cancer [11]. Furthermore, seven lncRNAs (MAGI2-AS3, GGTA1P, NAP1L2, CRABP2, SYNPO2, MKI67, and COL4A6) detected to be associated with TNBC prognosis, can be promising biomarkers [12]. Advancements in microarray and high throughput sequencing technologies have provided efficient tools to help in developing more reliable biomarkers for diagnosis, survival and prognosis [13,14]. However, the predictive power of a single gene biomarker may be insufficient. Emerging studies have found that gene signatures, including several genes, may be better alternatives [15]. To the best of our knowledge, the studies about multi-gene prognostic signatures in TNBC are very few, and the functions and mechanisms of mRNAs in TNBC remain to be further explored. Thus, it is necessary to identify more sensitive and efficient mRNA signatures for TNBC prognosis.
In this study, we first identified differentially expressed genes (DEGs), using 1109 BC samples and 113 matched non-cancerous samples from The Cancer Genome Atlas (TCGA). We identified ten hub genes associated with the cell cycle by functional enrichment analysis, protein-protein interaction (PPI) network and survival analysis. In addition, we developed a novel six-gene signature that could effectively predict TNBC survival.

Collection of clinical specimen data from the TCGA and GEO databases
The mRNA expression profiles and corresponding clinical information of breast cancer patients were downloaded from the Cancer Genome Atlas (TCGA) and gene expression omnibus (GEO) databases. We collected 1109 samples with gene expression data, containing 1109 BC tumour tissues samples and 113 normal tissue samples from TCGA database. After removing patients with incomplete information, we were left with117 TNBC samples and 970 non-TNBC samples. We collected 270 samples with 58 normal breast tissue samples and 212 TNBC tissue samples from the GEO dataset of the NCBI GEO database (GSE31519, GSE9574, GSE20194, GSE20271, GSE45255, and GSE15852).

Identification of differentially expressed genes
First, we merged the RNA-sequencing (RNA-seq) dataset files into a matrix file using the Perl language merge script. The gene name was converted from an Ensembl id to a gene symbol via the Ensembl database. Finally, the "edgeR" and "heatmap"R package were used to screen for differential genes between 117 TNBC and 970 other breast cancer patient subtypes and to map volcanoes. | log FC | > 1.0 and P < 0.05 were considered as the threshold value.

Functional and pathway enrichment analysis
Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of DEGs were performed using Database for Annotation, Visualization and Integration Discovery, DAVID version 6.8 [16]. P < 0.05 was chosen as the cutoff criterion. GO is a set of unified vocabulary to describe molecular functions (MF), biological processes (BP) and cellular components (CC) of biology, whereas KEGG analysis was performed to aid understanding of the signalling pathways involving DEGs.

PPI network construction and modules selection
A PPI network of differential genes was constructed, using STRING version 10.5 to evaluate information on protein-protein interactions [17]. Using the Molecular Complex Detection (MCODE) plug-in in Cytoscape 3.7.0, a visualization tool for integrating many molecular states such as expression level and interaction information into a unified conceptual framework [18], the PPI network module with densely connected regions was obtained (Degree cut-off > 15) [19].

Survival analysis
Clinical characteristic information for breast cancer was downloaded from TCGA. After removing samples with incomplete clinical overlapping DEG data, samples from 117 TNBC patients were used for further analysis. Univariate and multivariate Cox model analyses were used to identify candidate genes that were significantly associated with overall survival (OS). Based on the expression level and coefficient (β) of each gene, calculated by multivariate Cox proportional hazards regression analysis, a novel reliable prognostic gene signature was established. These TNBC patient samples were further divided into low-or high-risk groups based on the median risk score as the cut-off point. Kaplan-Meier curves were used to assess the prognostic value of the risk score. In addition, a time-dependent receiver operating characteristic (ROC) curve analysis, using the R package "survivalROC" was constructed to assess the predictive accuracy of the gene signature for time-dependent cancer death [20].
The area under the curve (AUC) was calculated to evaluate the predictive ability of the gene signature for clinical outcomes.

GO term and KEGG pathway enrichment analysis of DEGs
GO function and KEGG pathway enrichment analysis were performed using DAVID to expound the biological functions of 755 DEGs (Additional file 2: Table S2). The BP results indicated that DEGs were mainly significantly enriched in mitotic nuclear division, sister chromatid cohesion, cell division (Fig. 2a). MF analysis showed that DEGs were significantly enriched in microtubule motor, chemokine and structural molecule activities (Fig. 2b). CC analysis showed that the DEGs were mainly enriched in the extracellular region, chromosome centromeric region and kinetochore (Fig. 2c). In addition, the most enriched KEGG pathways were PPAR signalling, AMPK signalling and oocyte meiosis pathways (Fig. 2d).

A cell cycle related module selection by PPI Network analysis
Protein interactions among overlapping DEGs were predicted with STRING tools. A total of 148 nodes and 477 edges were displayed in the PPI network ( Fig. 3) with PPI-enrichment P value < 1.0e−16. A PPI network for DEG subsets with a combined score > 0.9 was constructed to determine the candidate hub genes. Based on the PPI network of the subsets, a module with an MCODE score of 42 and 45 nodes was identified ( Fig. 4a), and functional enrichment analyses showed that the genes in this module were mainly associated with the cell cycle and mitosis ( Fig. 4b and Table 1). BP analysis showed that these genes were significantly enriched in microtubule-based movement, mitotic sister chromatid segregation, mitotic metaphase plate congression, cell division, and mitotic cytokinesis. For CC analysis, these genes were significantly enriched in the condensed nuclear chromosome outer kinetochore, kinetochore, and spindle midzone. MF analysis showed the genes were significantly enriched in ATP binding, microtubule motor activity, single-stranded DNA binding, and DNA replication origin binding. In addition, the results of KEGG pathway enrichment analysis suggested that the pathways were enriched as follows: cell cycle, progesterone-mediated oocyte maturation, and oocyte meiosis. As a result, the eight genes correlated with cell cycle were selected as hub genes, which were CCNA2, CCNB2, CDC20, BUB1, TTK, CENPF, CENPA and CENPE (Table 2). Their expression levels were validated in 117 TNBC samples and 113 normal controls with breast cancer mRNA data from TCGA. As shown in Fig. 5, the eight mRNAs were significantly increased in TNBC compared with 113 normal control tissues (P < 0.001). We validated on the GEO database that the eight mRNA were also significantly increased compared with normal control tissues in TNBC (P < 0.001) (Additional file 3: Fig. S1). Using Cox proportional hazards regression model, we analysed the genes in the module, but no significant gene signature was established to predict overall survival.

Construction of a six-mRNA signature for survival prediction
A total of 16 out of 755 DEGs were significantly correlated with survival time (P < 0.05) and identified by the univariate Cox proportional hazards regression model (Additional file 2: Table S3). Additionally, a prognostic gene signature, composed of six genes, was developed after using the multivariate Cox proportional hazards regression model. The genes include transmembrane protein 252 (TMEM252), collagen type IX alpha 3 chain (COL9A3), proline rich protein BstNI subfamily 2 (PRB2), single-pass membrane protein with coiled-coil domains 1 (SMCO1), involucrin (IVL), and submaxillary gland androgen regulated protein 3B (SMR3B) ( Table 3). Patients were divided into low-and high-risk groups by the median risk score (1.070) (risk score = expression of SMR3B × 1.2141 + expression of TMEM252 × 1.6187 + expression of PRB2 × 1.4416 + expression of PRB2 × 2.0147 + expression of SMCO1 × 1.1471 + expression of COL9A3 × − 0.6101). The sixgene-based risk score distribution was presented in Fig. 6a. A highly significant difference in overall survival (OS) was detected between high-and low-risk groups (P < 0.0001) as shown in Fig. 6b. Moreover, the survival rate of the high-risk group was significantly much lower than for the low-risk group as depicted by Kaplan-Meier analysis in Fig. 6c (P < 0.0001). Time dependent ROC curve revealed that the prognostic signature presented a good performance in survival prediction, as shown in Fig. 6d and that the AUC was 0.929for 3 years OS and 0.902 for 5 years. Expression levels of the six genes in low-and high-risk groups are shown in Fig. 6e.

6-mRNA signature act as an independent prognostic indicator
Using univariate and multivariate Cox regression analyses, we investigated whether the prognostic values of the six mRNA were independent of clinicopathological factors. Univariate Cox regression model showed that the risk score, race, TNM stage, N status, M status, tumour status, and radiation were significantly related to the patients' overall survival in patients with TNBC (Table 4). In addition, multivariate Cox analysis indicated that the risk score and N stage still had remarkable independent prognostic values, with P = 0.005 and 0.025, respectively (Table 4). These results indicate that the 6-mRNA risk score was an independent prognostic indicator that can effectively predict the prognosis of TNBC patients.

Discussion
TNBC is characterized as a complex and aggressive disease with poor survival rates compared with other subtypes. Only 30% to 45% of TNBC patients achieve a complete pathological response and survival rates similar to other breast cancer subtypes [21]. The poor prognosis of patients diagnosed with TNBC is mainly due to a lack of effective targets for treatment. Therefore, there is an urgent need for more effective therapeutic targets to improve TNBC prognosis.
Misregulation of the cell cycle is a hallmark of cancer [22], disorders in mechanisms of cell cycle monitoring and proliferation cause tumour cell growth and tumour cell-specific phenomena. However, it remains unclear if misregulation of periodic mRNAs bears significance in TNBC patient pathogenesis. In this study, a total of 755 DEGs involved in TNBC were screened out from TCGA database, including 590 up-regulated and 165 downregulated genes. We then built related PPI networks of these DEGs and identified a significant module related to cell cycle, including several key DEGs in the regulatory network of TNBC patients. Subsequently, we identified eight periodic core genes (CCNA2, CCNB2, CDC20, BUB1, TTK, CENPF, CENPA, and CENPE) in the PPI network with higher capacity for PPIs. Coincidentally, all of them were up-regulated genes in TNBC (Fig. 5). CCNA2 (CyclinA2) and CCNB2 (CyclinB2) are members of the cyclin family of proteins that play key roles in the progression of G2/M transition, and have been reported to be the risk factors for resistance and recurrence [23][24][25]. Importantly, CCNA2, CCNB2, CDC20, BUB1, TTK, CENPA, and CENPE have been reported to be potential therapeutic targets for TNBC [26][27][28][29], and TTK inhibitors are currently being evaluated as anticancer   Clinical pathological features (Additional file 2: Table S4) are the proper prognostic references for TNBC patients. However, recent studies have demonstrated that clinical predictors are insufficient to precisely predict patient disease outcomes. The mRNA prognostic biomarker has the robust capacity of predicting the survival status of cancer patients. For example, Papadakis et al. [30] confirmed that mRNA BAG-1 acts as a biomarker in early breast cancer prognosis, Zheng et al. [31] found that CBX2 is a potential prognostic biomarker and therapeutic target for breast cancer.
However, it is insufficient as the single gene marker to independently predict patient survival. Because a single gene is easily affected by various factors, it is difficult to provide a stable and effective prediction effect. Therefore, we used Cox model analysis to construct a gene signature that includes several genes to enhance prognostic prediction efficiency and sensitivity to TNBC. It has been widely confirmed that combined genetic models are superior to previous single gene markers in disease prediction and diagnoses [32].
In this study, we constructed a six-mRNA (TMEM252, PRB2, SMCO1, IVL, SMR3B and COL9A3) signature for efficient and sensitive prognosis of TNBC patients. A previous study reported that COL9A3 potentially contributes to the pathogenesis of canine mammary tumours [33]. In another study, using RNA-seq to identify diabetic nephropathy, the expression of TMEM252, increased in diabetic patients relative to wild-type controls [34], but we have not found any relevant studies of TMEM252 in tumours. PRB2 is a key factor in regulating ER gene expression. In MCF-7 cells, PRB2 can interact with ER-beta to interfere with ER-beta shuttle between nuclear and cytoplasm [35], whereas ER-α gene inactivation is mediated by PRB2 in ER-negative breast cancer cells [36]. These findings suggest that PRB2 may be considered a promising target for TNBC therapy. Only one NCBI article was found to study the function of the single-pass membrane protein with coiled-coil domains 1 (SMCO1), which may contribute to hepatocyte proliferation and have the potential to promote liver repair and regeneration [37]. However, we have not found any research on SMCO1 in breast cancer; we speculate that it may also play an important role in breast cell proliferation. Additionally, we are not aware of any specific study on SMR3B in tumours, but SMR3B amplification has been detected in osteopontin (OPN)-positive hepatocellular carcinoma [38]. Involucrin (IVL), a component of keratinocyte crosslinked envelope, is found in the cytoplasm and crosslinked with membrane proteins by transglutaminase. This gene is mapped to 1q21, among calpactin I light chain, trichohyalin, profillaggrin, loricrin, and calcyclin. However, to our knowledge, there is no research on IVL in TNBC.
As far as we know, this is the first established 6-mRNA signature for the prediction of OS time in TNBC, and we have demonstrated the independent prognostic value of this 6-mRNA signature in TNBC.

Conclusions
In summary, through bioinformatic analysis, we identified eight hub genes, correlated with cell cycle, that might be tightly correlated with TNBC pathogenesis. Besides, we constructed a 6-mRNA signature which may act as a potential prognostic biomarker in patients with TNBC, and the prognostic model presented a good performance in OS prediction at 3 and 5 years. These findings will provide some guidance for future TNBC prognosis

Additional files
Additional file 1: Table S1. Differential expressed mRNAs between 1109 breast cancer tissue samples and 113 normal tissue samples. Differential expressed mRNAs between 117 TNBC and 970 non-TNBC breast cancer samples. 590 up-regulated overlapped DEGs and 165 down-regulated overlapped DEGs.
Additional file 2: Table S2. Functional and pathway enrichment analysis of DEGs in TNBC.