Abstract Background Triple-negative breast cancer (TNBC) is a subtype of breast cancer with poor clinical outcome and limited treatment options. Lacking molecular targets, chemotherapy is the main adjuvant treatment for TNBC patients. Materials and methods To explore potential therapeutic targets for TNBC, we analyzed three microarray datasets ([37]GSE38959, [38]GSE45827, and [39]GSE65194) derived from the Gene Expression Omnibus (GEO) database. The GEO2R tool was used to screen out differentially expressed genes (DEGs) between TNBC and normal tissue. Gene Ontology function and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis were performed using the Database for Annotation, Visualization and Integrated Discovery to identify the pathways and functional annotation of DEGs. Protein–protein interaction of these DEGs was analyzed based on the Search Tool for the Retrieval of Interacting Genes database and visualized by Cytoscape software. In addition, we used the online Kaplan–Meier plotter survival analysis tool to evaluate the prognostic value of hub genes expression in breast cancer patients. Results A total of 278 upregulated DEGs and 173 downregulated DEGs were identified. Among them, ten hub genes with a high degree of connectivity were picked out. Overexpression of these hub genes was associated with unfavorable prognosis of breast cancer, especially, CCNB1 overexpression was observed and indicated poor outcome of TNBC. Conclusion Our study suggests that CCNB1 was overexpressed in TNBC compared with normal breast tissue, and overexpression of CCNB1 was an unfavorable prognostic factor of TNBC patients. Further study is needed to explore the value of CCNB1 in the treatment of TNBC. Keywords: triple-negative breast cancer, hub genes, expression profiling data, CCNB1 Introduction Triple-negative breast cancer (TNBC) is defined as a subtype of breast cancer which lacks expression of estrogen receptor (ER) and progesterone receptor (PR) and demonstrates no amplification of human epidermal growth factor receptor 2 (HER2). This subset accounts for ~12%–17% of all invasive breast cancers.[40]^1 TNBC is more frequently diagnosed in younger women and behaves more aggressively in clinical behaviors. Patients with TNBC are more likely to develop relapse and visceral metastasis than other subtypes of breast cancer.[41]^2^–[42]^5 Lacking molecular targets, patients diagnosed with TNBC cannot be treated with endocrine therapy or HER2-targeted therapy. Chemotherapy is currently the main adjuvant treatment for TNBC patients.[43]^1 Unfortunately, many tumors are resistant to chemotherapy and relapse or metastasize quickly after adjuvant treatment.[44]^6^,[45]^7 Up to date, TNBC is still a disease with poor outcome and limited treatment options. Hence, it is urgent and necessary to explore novel therapeutic targets for TNBC. In this study, we tried to detect novel indicators of poor prognosis in TNBC patients and endeavor to provide potential therapeutic targets for this challenging disease. To detect the differentially expressed genes (DEGs) between TNBC and healthy human breast tissue, bioinformatics methods were used to analyze the gene expression profiling data downloaded from the Gene Expression Omnibus (GEO) database. Gene Ontology (GO) functional annotation analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed for the screened DEGs. Then, we established a protein–protein interaction (PPI) network to identify hub genes related to TNBC. The survival analysis of these hub genes was performed using the online database Kaplan–Meier plotter. Materials and methods Data source The gene expression datasets analyzed in this study were obtained from the GEO database ([46]https://www.ncbi.nlm.nih.gov/geo/). A total of 1,821 series about human breast cancer were retrieved from the database. After a careful review, three gene expression profiles ([47]GSE38959, [48]GSE45827, and [49]GSE65194) were selected. Among them, [50]GSE38959 was based on the Agilent [51]GPL4133 platform (Agilent-014850 Whole Human Genome Microarray 4×44K G4112F), and [52]GSE45827 and [53]GSE65194 were based on platform [54]GPL570 ([HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array). All of the data were freely available online, and this study did not involve any experiment on humans or animals performed by any of the authors. Data processing of DEGs The GEO2R online analysis tool ([55]https://www.ncbi.nlm.nih.gov/geo/geo2r/) was used to detect the DEGs between TNBC and normal samples, and the adjusted P-value and |logFC| were calculated. Genes that met the cutoff criteria, adjusted P<0.05 and |logFC|≥2.0, were considered as DEGs. Statistical analysis was carried out for each dataset, and the intersecting part was identified using the Venn diagram webtool ([56]bioinformatics.psb.ugent.be/webtools/Venn/). GO and KEGG pathway analysis of DEGs GO analysis is a common useful method for large scale functional enrichment research; gene functions can be classified into biological process (BP), molecular function (MF), and cellular component (CC). KEGG is a widely used database which stores a lot of data about genomes, biological pathways, diseases, chemical substances, and drugs. GO annotation analysis and KEGG pathway enrichment analysis of DEGs in this study was performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID) tools ([57]https://david.ncifcrf.gov/). P<0.01 and gene counts≥10 were considered statistically significant. PPI network construction and hub gene identification The Search Tool for the Retrieval of Interacting Genes (STRING) database ([58]http://string-db.org/) is designed to analyze the PPI information. To evaluate the potential PPI relationship, the DEGs identified previously were mapped to the STRING database. The PPI pairs were extracted with a combined score>0.4. Subsequently, the PPI network was visualized by Cytoscape software ([59]www.cytoscape.org/). Nodes with higher degree of connectivity tend to be more essential in maintaining the stability of the entire network. CytoHubba, a plugin in cytoscape, was used to calculate the degree of each protein node. In our study, the top ten genes were identified as hub genes. Survival analysis of hub genes The Kaplan–Meier plotter ([60]http://kmplot.com/analysis/) is an online tool applied to assess the effect of 54,675 genes on survival using 10,461 cancer samples (5,143 breast, 1,816 ovarian, 2,437 lung, and 1,065 gastric cancer). The Kaplan–Meier plotter mRNA breast cancer database was applied to evaluate the prognostic values of hub genes in breast cancer patients, especially in TNBC patients. In our study, TNBC patients were screened out based on ER, PR, and HER-2 negative expression. Probes of genes were selected based on the “only JetSet best probe set,” and the desired probe IDs for each gene are shown in [61]Table S1. For each gene, cancer patients were divided into two groups according to the median values of mRNA expression. P<0.01 was considered to indicate a statistically significant result. Results Identification of DEGs Three gene expression profiles ([62]GSE38959, [63]GSE45827, and [64]GSE65194) were selected in this study. Among them, [65]GSE38959 contained 30 TNBC samples and 13 normal samples, and [66]GSE45827 and [67]GSE65194 included 41 TNBC specimens and eleven normal breast specimens respectively ([68]Table 1). Based on the criteria of P<0.05 and |logFC|≥2, a total of 852 DEGs were identified from [69]GSE38959, including 515 upregulated genes and 337 downregulated genes. In gene chip [70]GSE45827, 2,995 DEGs were identified; 2,117 genes were upregulated, and 878 genes were downregulated. And from [71]GSE65194, 3,031 DEGs including 2,130 upregulated genes and 901 downregulated genes were identified. All DEGs were identified by comparing TNBC samples with normal breast samples. Subsequently, Venn analysis was performed to get the intersection of the DEG profiles ([72]Figure 1). Finally, 451 DEGs were significantly differentially expressed among all three groups, of which 278 were significantly upregulated genes and 173 were downregulated. Table 1. Statistics of the three microarray databases derived from the GEO database Dataset ID TNBC Normal Total number [73]GSE38959 30 13 43 [74]GSE45827 41 11 52 [75]GSE65194 41 11 52 [76]Open in a new tab Abbreviations: GEO, Gene Expression Omnibus; TNBC, triple-negative breast cancer. Figure 1. [77]Figure 1 [78]Open in a new tab Venn diagram of DEGs common to all three GEO datasets. Notes: (A) Downregulated genes. (B) Upregulated genes. Abbreviations: DEG, differentially expressed gene; GEO, Gene Expression Omnibus. Functional enrichment analyses of DEGs GO function and KEGG pathway enrichment analysis for DEGs were performed using the DAVID ([79]Table 2). The enriched GO terms were divided into CC, BP, and MF ontologies. The results of GO analysis indicated that DEGs were mainly enriched in BPs, including sister chromatid cohesion, microtubule-based movement, anaphase-promoting complex-dependent catabolic process, and extracellular matrix (ECM) organization. MF analysis showed that the DEGs were significantly enriched in microtubule binding, transcriptional activator activity, RNA polymerase II core promoter proximal region sequence-specific binding, ATPase activity, and microtubule motor activity. For the cell component, the DEGs were enriched in condensed chromosome kinetochore, microtubule, kinetochore, and spindle. In addition, the results of KEGG pathway analysis showed that DEGs were mainly enriched in pathways in cancer, small cell lung cancer, and ECM–receptor interaction. Table 2. Significantly enriched GO terms and KEGG pathways of DEGs Category Term Description Count P-value BP term GO:0007062 Sister chromatid cohesion 23 4.30E–15 BP term GO:0007018 Microtubule-based movement 11 2.59E–05 BP term GO:0031145 Anaphase-promoting complex-dependent catabolic process 10 1.24E–04 BP term GO:0030198 ECM organization 14 1.03E–03 CC term GO:0000777 Condensed chromosome kinetochore 23 6.49E–17 CC term GO:0000776 Kinetochore 17 7.09E–11 CC term GO:0005819 Spindle 14 6.63E–06 CC term GO:0005874 Microtubule 19 5.31E–04 MF term GO:0008017 Microtubule binding 20 7.71E–07 MF term GO:0003777 Microtubule motor activity 11 2.27E–05 MF term GO:0001077 Transcriptional activator activity, RNA polymerase II core promoter proximal region sequence-specific binding 15 1.90E–03 MF term GO:0016887 ATPase activity 12 5.11E–03 KEGG pathway hsa05222 Small cell lung cancer 11 1.76E–04 KEGG pathway hsa05200 Pathways in cancer 26 1.95E–04 KEGG pathway hsa04512 ECM–receptor interaction 11 2.14E–04 [80]Open in a new tab Abbreviations: BP, biological process; CC, cellular component; DEG, differentially expressed gene; ECM, extracellular matrix; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function. PPI network construction and hub gene identification Protein interactions among the DEGs were predicted with STRING tools. A total of 111 nodes and 1,365 edges were involved in the PPI network, as presented in [81]Figure 2. The top ten genes evaluated by connectivity degree in the PPI network were identified ([82]Table 3). The results showed that cyclin-dependent kinases 1 (CDK1) was the most outstanding gene with connectivity degree=64, followed by cyclin B1 (CCNB1; degree=61), baculoviral IAP repeat containing 5 (BIRC5; degree=60), aurora kinase A (AURKA; degree=58), polo-like kinase 1 (PLK1; degree=56), mitotic arrest deficient 2-like 1 (MAD2L1; degree=54), BUB1 mitotic checkpoint serine/threonine kinase B (BUB1B; degree=54), nuclear division cycle 80 (NDC80; degree=53), budding uninhibited by benzimidazoles 1 (BUB1; degree=52), and kinesin family member 11 (KIF11; degree=52). All of these hub genes were upregulated in TNBC. Figure 2. [83]Figure 2 [84]Open in a new tab Protein–protein interaction network constructed with the differentially expressed genes. Note: Red nodes represent upregulated genes, and blue nodes represent downregulated genes. Table 3. Top ten hub genes with higher degree of connectivity Gene symbol Gene description Degree CDK1 Cyclin-dependent kinases 1 64 CCNB1 Cyclin B1 61 BIRC5 Baculoviral IAP repeat containing 5 60 AURKA Aurora kinase A 58 PLK1 Polo-like kinase 1 56 MAD2L1 Mitotic arrest deficient 2-like 1 54 BUB1B BUB1 mitotic checkpoint serine/threonine kinase B 54 NDC80 Nuclear division cycle 80 53 BUB1 Budding uninhibited by benzimidazoles 1 52 KIF11 Kinesin family member 11 52 [85]Open in a new tab Survival analysis of ten hub genes To investigate the prognostic values of the ten potential hub genes, the Kaplan–Meier plotter bioinformatics analysis platform was used. A total of 1,402 breast cancer patients were available for the analysis of overall survival. We found that high expression of these hub genes was associated with unfavorable overall survival of breast cancer patients ([86]Figure 3). Figure 3. [87]Figure 3 [88]Open in a new tab Kaplan–Meier overall survival analyses for the top ten hub genes expressed in breast cancer patients. Note: See [89]Table 3 for gene description. However, only overexpression of CCNB1 was an unfavorable prognostic factor of relapse-free survival in TNBC patients (HR=2.12; 95% CI: 1.2–3.72; P=0.0078; n=255). There were not enough incidents for overall survival analysis ([90]Figure 4). Figure 4. Figure 4 [91]Open in a new tab Kaplan–Meier relapse-free survival analyses for CCNB1 expression in TNBC patients. Abbreviations: CCNB1 cyclin B1; TNBC, triple-negative breast cancer. Discussion Breast cancer is a heterogeneous disease, and the histopathological features and clinical behaviors are distinct among subtypes. TNBC is a unique subtype of breast cancer with poor prognosis. Patients with TNBC have an increased likelihood of relapse and visceral metastasis. Due to lacking a therapeutic target, patients with TNBC could not benefit from endocrine therapy or HER2-targeted therapy, and chemotherapy is currently the mainstay of adjuvant treatment. However, TNBC patients are more likely to develop chemoresistance. Hence, it is crucial to identify new specific targeted therapies for TNBC. In the present study, gene expression and protein–protein expression analysis based on publicly available databases was performed to identify potential key genes correlated with TNBC. DEGs between TNBC and healthy human breast tissue were screened out based on gene expression profiling data from the GEO database. Totally, we identified 278 upregulated DEGs and 173 downregulated DEGs. These DEGs were associated with the GO BP terms such as condensed chromosome kinetochore, sister chromatid cohesion, kinetochore, and microtubule binding, and significantly enriched in the KEGG terms small cell lung cancer, pathways in cancer, and ECM–receptor interaction. A PPI network was constructed to investigate the interrelationship of the DEGs, and ten hub genes were identified, including AURKA, BIRC5, BUB1B, BUB1, CCNB1, CDK1, KIF11, MAD2L1, NDC80, and PLK1. All of these genes were upregulated in TNBC. Finally, the Kaplan–Meier plotter online tool was applied to predict the relationship between the expression of hub genes and prognosis of TNBC patients. Based on the Kaplan–Meier plotter, overexpression of all the above genes was related to unfavorable prognosis of breast cancer patients. However, only overexpression of CCNB1 was an unfavorable prognostic factor of TNBC patients. CCNB1, also known as cyclin B1, is a key modulator in controlling cell proliferation.[92]^8 Some research has demonstrated that CCNB1 is involved in apoptosis, chemoresistance, and epithelial mesenchymal transitions of tumor cells.[93]^9^,[94]^10 Overexpression of cyclin B1 has been reported in many tumors, such as colorectal cancer, gastric cancer,[95]^11 pancreatic carcinoma,[96]^12 and lung carcinoma.[97]^13 Some of these studies suggested that the overexpression of cyclin B1 may be associated with the poor prognosis of these malignant diseases. For breast cancer, a lot of studies have shown that cyclin B1 overexpression was associated with aggressive clinical behaviors and was an independent prognostic factor. Aaltonen et al[98]^14 showed that cyclin B1 overexpression was correlated with an aggressive phenotype and was significantly associated with shorter overall survival and metastasis-free survival in breast cancer patients. Ding et al[99]^15 reported that a high level of CCNB1 was closely associated with hormone therapy resistance and poor recurrence-free survival, disease-free survival, and distant metastasis-free survival of ER+ breast cancer patients. And a meta-analysis by Sun et al[100]^16 suggested that cyclin B1 overexpression might be an independent potential prognostic marker for disease-specific survival and disease-free survival of breast cancer. TNBC are usually high grade tumors with primitive features, suggesting that cyclin B1 may overexpress in TNBC. Agarwal et al[101]^17 reported that cyclin B1 was expressed at a significantly higher level in TNBC cell lines than other subtypes. In our study, cyclin B1 was overexpressed in TNBC compared to normal breast tissue, and overexpression of cyclin B1 was correlated with unfavorable relapse-free survival of TNBC patients. Therefore, cyclin B1 may be a prognostic factor and potential therapeutic target for TNBCs. Except for CCNB1, we detected other nine hub genes associated with breast cancer, including CDK1, AURKA, BIRC5, MAD2L1, BUB1B, BUB1, PLK1, KIF11, and NDC80. Most of them were reported as an essential factor involved in cell division and proliferation. Proteins encoded by AURKA, BUB1, BUB1B, PLK1, and CDK1 are all serine/threonine kinases involved in the regulation of the cell cycle,[102]^18 and overexpression of these genes has been detected in various human cancers and correlated with their prognosis. Roylance et al[103]^19 reported that a high AURKA expression level was significantly associated with poorer clinical outcome in breast cancer patients. In Sotiriou et al’s study,[104]^20 BUB1 was upregulated and correlated with a poor clinical prognosis in breast cancer patients. Many studies have shown an association between PLK1 overexpression and poor clinical prognosis, and suggested that inhibition of PLK1 may be a potential therapy for cancer treatment.[105]^21^,[106]^22 For CDK1, many research studies have reported its overexpression in cancers and that it acts as an adverse prognostic factor, and many kinds of CDK inhibitors have been developed.[107]^23 In our study, BIRC5, KIF11, MAD2L1, and NDC80 were overexpressed in breast cancer compared to normal breast tissues, and overexpression of these genes was significantly correlated with unfavorable clinical outcome in breast cancer patients. The results of our research were consistent with other studies.[108]^24^–[109]^27 However, the role of these genes in TNBC is not clear and further study is needed. Conclusion Our bioinformatics analysis identified 451 DEGs between TNBCs and normal breast tissues based on the gene expression datasets obtained from the GEO database. Among them, ten hub genes might be the core genes of breast cancer, including AURKA, BIRC5, BUB1B, BUB1, CCNB1, CDK1, KIF11, MAD2L1, NDC80, and PLK1. All of them were upregulated in breast cancer, and overexpression of these genes was associated with unfavorable clinical outcome in breast cancer patients. In TNBC patients, CCNB1 overexpression is an unfavorable prognostic factor. Further study is needed to confirm the results of our research. Anyway, CCNB1 may be a potential target for TNBC therapy. Supplementary material Table S1. The desired probes of hub genes in the Kaplan–Meier plotter database Gene symbol Probe ID CDK1 203213_at CCNB1 228729_at BIRC5 202094_at AURKA 208079_s_at PLK1 202240_at MAD2L1 203362_s_at BUB1B 203755_at NDC80 204162_at BUB1 209642_at KIF11 204444_at [110]Open in a new tab Acknowledgments