Identification of hub genes and small molecule therapeutic drugs related to breast cancer with comprehensive bioinformatics analysis

Download of datasets and identification of DEGs

Five independent gene expression profiles (GSE3744, GSE21422, GSE42568, GSE61304, and GSE65194) based on GPL570 Platform (Affymetrix Human Genome U133 Plus 2.0 Array) were downloaded from GEO database to identify DEGs. After normalization between arrays, we investigated DEGs among each dataset with the threshold of |log2FoldChange (log2FC)|> 1 and adj.P.Val < 0.05 using LIMMA package of R (Ritchie et al., 2015).

For validation, GEPIA2 online tool (http://gepia2.cancer-pku.cn/#index) was used to analyze the differential expression of TCGA Breast invasive carcinoma (BRCA) RNA-seq dataset composed of 1,085 tumor samples and 112 normal samples according to the cut-off standard (|log2FC|>1 and q-value <0.05). Lastly, the common results of TCGA BRCA dataset and GEO datasets were selected as the overlapping DEGs of breast cancer, which could reduce the influences resulted from the heterogeneity of the different datasets. Venn diagram and volcano plot were drawn by VennDiagram and other packages of R. Table 1 listed the details of datasets.

Protein–protein interaction (PPI) network integration and hub genes screening

STRING (version 11.0; https://string-db.org/) is a biological database designed to analyze functional interactions between proteins (Szklarczyk et al., 2019). In this study, we used STRING to construct a PPI network with overlapping DEGs under the premise of an Interaction score of 0.7. Then, we utilized the cytoHubba plug-in of Cytoscape (version 3.7.2) which provided the calculated results by maximal clique centrality (MCC), maximum neighborhood component (MNC) and Degree methods to identify hub genes from the PPI network (Chin et al., 2014).

Verification of hub genes

The Oncomine database (https://www.oncomine.org) was used to verify the mRNA expression of hub genes with the threshold of P <0.05 and fold change >2. Next, the Human Protein Atlas database (HPA, https://www.proteinatlas.org) was used to validate the protein expression of genes by immunohistochemistry data.

BC-GenExMiner (http://bcgenex.centregauducheau.fr/BC-GEM/GEM-Accueil.php) is a statistical mining tool that contains published breast cancer transcription data (10,716 DNA microarray samples and 4,712 RNA-seq samples) (Jézéquel et al., 2012). We performed correlation analysis between hub genes in breast cancer using BC-GenExMiner and GEPIA2 online tool.

GO functional and KEGG pathway enrichment analysis

ClusterProfiler package of R can automatically classify biological terms and gene clusters (Yu et al., 2012). To elucidate the biological characteristics of breast cancer-related genes, we performed GO functional and KEGG pathway enrichment analysis by ClusterProfiler with p-value <0.05 and q-value <0.05.

Analysis of genetic alterations of hub genes

The cBioPortal for Cancer Genomics (cBioPortal; http://cbioportal.org) provides online resources for the exploration, visualization and analysis of multidimensional cancer genomics data (Gao et al., 2013). In this study, 6618 breast cancer samples from 13 related reports in cBioPortal were used as research materials to explore genetic alterations connected with the selected hub genes. Afterward, we utilized COSMIC (https://cancer.sanger.ac.uk/cosmic), the most comprehensive resource for studying somatic mutation information in human cancer, to analyze hub genes alterations in breast cancer (Forbes et al., 2010).

Evaluation of clinicopathological characteristics and immune infiltration

We used BC-GenExMiner to analyze the correlations between hub genes expression and clinicopathological variables such as Oestrogen receptor status (ER), Progesterone receptor status (PR), HER2 receptor status (HER2), Nodal status (N), Scarff Bloom & Richardson grade status (SBR), Nottingham Prognostic Index status (NPI), Age status, P53 status, Basal-like and Triple negative breast cancer (TNBC) subtypes. P < 0.05 was considered to be statistically significant.

TIMER (https://cistrome.shinyapps.io/timer/) is a comprehensive resource for systematic analysis of tumor-infiltrating immune cells across 32 different cancers from TCGA database (Li et al., 2017). In this experiment, we estimated the associations between hub genes expression and immune cell populations (B Cell, CD8+ T Cell, CD4+ T Cell, Macrophage, Neutrophil, and Dendritic Cell) in breast cancer using TIMER.

Prediction of TFs-hub genes and miRNAs-hub genes interaction

TRRUST (version 2, https://www.grnpedia.org/trrust/), a manually curated database of human and mouse transcriptional regulatory networks, was used to explore the TFs of hub genes (Han et al., 2018). Meanwhile, we used Encyclopedia of RNA Interactomes Platform (ENCORI, http://starbase.sysu.edu.cn/) to unearth the miRNAs targeting hub genes (Li et al., 2014). In addition, we utilized GEO2R online tool (https://www.ncbi.nlm.nih.gov/geo/geo2r/) to analyze breast cancer differentially expressed miRNAs (DEmiRNAs) in GSE97811 dataset including 45 tumor tissues and 16 normal tissues. Then, we got miRNAs-hub genes interaction combined with miRNAs-target genes analysis and DEmiRNAs data. Finally, Cytoscape was used to visualize transcriptional and post-transcriptional regulatory networks.

Survival analysis

In this section, we applied the Kaplan–Meier Plotter (http://kmplot.com/analysis/) to evaluate prognostic information of previously identified hub genes and important reporter regulatory molecules (Nagy et al., 2018). The expression values of these genes were split into either high (expression value ≥ median) or low (expression value <median). Hazard ratio (HR) was calculated to evaluate the association between genes expression and survival, and p < 0.05 was considered statistically significant.

Small molecule drugs analysis

The DGIdb online tool (http://www.dgidb.org/)—an available resource containing drug-gene interaction information from more than 30 databases—was used to screen antineoplastic drugs targeting hub genes (Cotto et al., 2018). We also downloaded the GSE85871 dataset, which is a gene expression data of MCF7 cells treated with 102 traditional Chinese medicine (TCM) ingredients recorded in the GEO database. LIMMA package of R was used to analyze the differential expression genes in each TCM ingredient treatment group compared with the untreated group (adj .p. val <  0.05). Then, we evaluated the reversal effects of each TCM ingredient on overlapping DEGs induced by breast cancer.

Dataset processing and DEGs acquisition

After data normalization, the black lines of gene expression box plots of all samples in each single dataset were almost on the same level, which was an important marker to predict the reliability and accuracy of the experimental results (Fig. S1). Next, DEGs (1,738 in GSE3744, 2,430 in GSE21422, 3,116 in GSE42568, 990 in GSE61304 and 4,181 in GSE65194) were identified (Figs. 2A–2E; Tables S1–S5), and 302 integrated DEGs (110 up- and 192 down-regulated) from 5 GEO datasets were found (Figs. 2F and 2G). Similarly, we identified 3,559 DEGs from TCGA BRCA dataset with the cut-off criteria of |Log2FC| > 1 and q-value < 0.05 (Table S6). Then, 283 overlapping DEGs which might play promoting or inhibitory roles in breast cancer progression were confirmed from the analysis results of GEO datasets and TCGA BRCA dataset (Figs. 2H and 2I).

PPI network construction and hub genes filtering

The PPI network around proteins encoded by 283 overlapping DEGs was constructed using STRING (Fig. S2A, Table S7). We found that 165 of the 283 overlapping DEGs were related to each other and were visualized using Cytoscape –165 nodes and 1,861 edges were included in PPI network and six hub genes (RRM2, CDC20, CCNB2, BUB1B, CDK1, and CCNA2) based on MCC, MNC and Degree methods were identified (Fig. S2B). Notably, these hub genes were all up-regulated in overlapping DEGs (Tables S1–S6) and might play important roles in the pathogenesis of breast cancer.

Verification of hub genes

Based on the large-scale breast cancer-related data in the Oncomine database, we confirmed that hub genes were significantly up-regulated in multiple cancer types, including Breast Cancer, Brain and CNS Cancer, Lymphoma, Lung Cancer, and so on (Fig. 3A). Also, immunohistochemistry staining data obtained from the HPA database demonstrated the up-regulated expression of proteins encoded by RRM2, CDC20, CCNB2, CDK1 and CCNA2 (Figs. 3B–3K). However, we did not find the association between BUB1B and breast cancer in HPA database. According to the current analysis, we predicted that BUB1B might also be associated with breast cancer, but experimental data were needed to confirm this specific connection. Meanwhile, CDK1, CCNA2, CCNB2, BUB1B and CDC20 were classified as cancer-related genes, and RRM2 was an FDA approved drug target. The data of BC-GenExMiner and GEPIA2 both confirmed a powerful correlation among hub genes (Fig. S3), suggesting that these genes might be the functional partners in breast cancer.

GO function enrichment and KEGG pathway analysis

The GO function annotations of overlapping DEGs were mainly classified into biological processes (BP), cell component (CC) and molecular function (MF). As for BP, up-regulated overlapping DEGs were significantly related to mitotic nuclear division, organelle fission, and regulation of chromosome segregation, which was consistent with the biological characteristics of the abnormally rapid proliferation of breast cancer cells. And the down-regulated genes were closely related to regulation of cellular response to growth factor stimulus, temperature homeostasis, retinoid metabolic process and regulation of vasculature development. Within CC, the up-regulated genes were remarkably correlated to spindle, chromosome, centromeric region, condensed chromosome kinetochore, and midbody, whereas the down-regulated genes were related to collagen-containing extracellular matrix and sarcolemma. MF analysis displayed the up-regulated genes were involved in microtubule binding, microtubule motor activity, and cyclin-dependent protein serine/threonine kinase regulator activity, whereas the down-regulated genes were mainly enriched in glycosaminoglycan binding, heparin binding, extracellular matrix structural constituent and growth factor binding (Figs. 4A and 4B; Tables S8 and S9).

KEGG pathway analysis showed the up- and down-regulated overlapping DEGs were all significantly attached to Cell cycle, Oocyte meiosis, Tyrosine metabolism, ECM-receptor interaction, Progesterone-mediated oocyte maturation, p53 signaling pathway, PPAR signaling pathway and Phenylalanine metabolism (Fig. 4C). Furthermore, the related pathways of hub genes included Cell cycle, Oocyte meiosis, Progesterone-mediated oocyte maturation, and p53 signaling pathway (Fig. 4C). Table S10 presented the detailed results of KEGG enrichment analysis.

Genetic alterations of hub genes

As a result, there were nearly 2.1% (CDC20, CDK1), 1.2% (RRM2), 0.9% (BUB1B), 0.8% (CCNA2), 0.7% (CCNB2) of breast cancer samples included in cBioPortal had genetic changes. Amplification was the most common genetic alterations among six hub genes. And deep deletion was another major genetic change among five hub genes (RRM2, CCNB2, BUB1B, CDK1, CCNA2). While genetic alterations in RRM2, CDC20, BUB1B, CDK1 and CCNA2 were related to missense mutation (Fig. S4A). In the query results of COSMIC database, we found that all hub genes had missense mutations. RRM2, CCNB2, CDK1, BUB1B and CDC20 experienced synonymous substitutions. CCNA2 and BUB1B also had frameshift deletions. Specifically, these hub genes had mutations such as A > G, A > T, C > T, C > G, G > A, G > T at nucleotide level (Fig. S4B).

Clinicopathological characteristics and immune infiltration

We investigated the relevance of six hub genes and clinicopathological features using BC-GenExMiner (Figs. S5–S7). Data analysis showed that higher expression of hub genes was found in higher NPI and SBR grade (p < 0.001). And the expression of these hub genes was significantly higher in ER-, PR-, HER2+, Nodal+, P53-mutated, Basal-like and TNBC clinical subtypes of breast cancer. Surprisingly, significantly increased expression of hub genes was found in patients not more than 51 years old (p < 0.001). Moreover, our current results demonstrated that hub genes were correlated with 6 types of immune cell infiltrates (B Cell, CD8+ T Cell, CD4+ T Cell, Macrophage, Neutrophil and Dendritic Cell) with various degrees based on the TIMER database (Fig. S8).

Prediction of TF-hub genes and miRNAs-hub genes interaction

We screened the potential regulatory relationships between TFs and hub genes via TRRUST database to further study the functional roles of hub genes. A total of 19 associations between 17 TFs and six hub genes were shown (Fig. 5A). E2F1 could activate the transcriptional process of RRM2 and CDK1. In contrast, TP53 inhibited the expression of CDK1 and CCNB2. Noticeably, PTTG1 gene, the transcriptional promoter of CDK1, was up-regulated in breast cancer; ZBTB16 gene was the transcriptional suppressor of CCNA2, and its expression was down-regulated in overlapping DEGs. Additionally, we obtained 273 targeted miRNAs with regulatory effects on hub genes using ENCORI online tool, and identified 103 DEmiRNAs (58 up- and 45 down-regulated; Table S11) from GSE97811 dataset with |log2FC|>1 and adj.P.Val < 0.05. As a result, 41 associations between 30 miRNAs and 4 hub genes (RRM2, CDK1, CCNA2, and BUB1B) were found from targeted miRNAs analysis and DEmiRNAs data (Fig. 5B). The top 2 hub genes with the most miRNAs targets were RRM2 and CCNA2. In addition, hsa-miR-340-5p, hsa-miR-130a-3p, hsa-miR-200b-3p, hsa-miR-200c-3p, hsa-miR-204-5p, hsa-miR-219a-5p, hsa-miR-27a-3p, hsa-miR-27b-3p, hsa-miR-301a-3p and hsa-miR-429 were the top 10 miRNAs with the most target genes. Unfortunately, we only used TRRUST and ENCORI databases, as well as limited samples of GSE97811 dataset to analyze TFs and miRNAs targeting hub genes, which may potentially limit the completeness of this study.

Survival analysis of hub genes

We then analyzed the prognostic information of six hub genes using Kaplan–Meier Plotter. The result demonstrated that breast cancer patients with higher hub genes expression had worse overall survival (OS), relapse-free survival (RFS) and distant metastasis-free survival (DMFS) (Figs. 6A–6R).

Furthermore, we performed prognostic analysis on the important reporter regulatory molecules and found that higher E2F1 and PTTG1 expression predicted worse OS. In contrast, higher expression of ZBTB16, hsa-miR-130a-3p and hsa-miR-204-5p was significantly associated with better OS in breast cancer (Fig. S9A). The remaining reporter regulatory molecules (e.g., TP53, hsa-miR-340-5p, hsa-miR-200b-3p, hsa-miR-200c-3p, hsa-miR-219a-5p, hsa-miR-27a-3p, hsa-miR-27b-3p, hsa-miR-301a-3p, hsa-miR-429) had no statistically significant correlation with OS (P >0.05; Fig. S9B).

Small molecule drugs analysis

Regarding six hub genes as potential therapeutic targets for breast cancer, we identified several antineoplastic drugs based on the DGIdb database. At present, only RRM2, CDK1 and CCNA2 were identified as tumor therapeutic targets. Therefore, we speculated the other three genes (CCNB2, BUB1B, CDC20) might be the novel targets in the future. Statistical analysis revealed that 21 candidate drugs such as Cladribine, Gallium nitrate, Dinaciclib, Alvocidib and Suramin targeted RRM2, CDK1 and CCNA2 (Table 2). More experimental data are needed to further confirm the potential of these drug candidates in the treatment of breast cancer.

In addition, we identified 8 TCM ingredients (1β-hydroxyalantolactone, Andrographolide, Berberine hydrochloride, Britanin, Hyodeoxycholic acid, Japonicone A, Nitidine chloride and Tanshinone IIA) that reversed breast cancer-induced gene expression from GSE85871 dataset. Japonicone A reversed the expression of 87 overlapping DEGs, including six hub genes, and its potential therapeutic effects on breast cancer were related to Cell cycle, Oocyte meiosis, p53 signaling pathway, Progesterone-mediated oocyte maturation, Viral carcinogenesis, HTLV-I infection, Pathways in cancer and TGF-beta signaling pathway (Fig. 7A). Also, Nitidine chloride, Berberine hydrochloride, 1β-hydroxyalantolactone, Britanin and Tanshinone IIA reversed the expression of 37, 35, 31, 30 and 19 overlapping DEGs, respectively, and the potential therapeutic effects of these ingredients on breast cancer were related to biological pathways such as Cell cycle, Oocyte meiosis, p53 signaling pathway and Progesterone-mediated oocyte maturation (Figs. 7B–7F). Details of TCM ingredients were shown in Table S12.