Human breast cancer microarray datasets were searched and downloaded from the NCBI-GEO database in March 2016. The key words ‘breast cancer’, ‘breast adenocarcinoma’ and ‘breast tumor’ were used to perform a specific search. The selection criteria for the datasets were as follows: i) All datasets were genome-wide; ii) the samples of each dataset included breast cancer patients and controls; iii) the samples in tumor and control group were from breast tissue; iv) the dataset included different stages of breast cancer; and v) raw data were available. Datasets were excluded if: i) The number of samples was <3 for cases or controls; and ii) severe RNA degradation or an insufficient number of detected genes. Based on the aforementioned criteria, five datasets were finally selected for the integrated analysis [GSE10810 (19), GSE16391 (20), GSE29431 (21), GSE42568 (22) and GSE61304 (23)]. The integrated datasets included 257 breast cancer patients and 98 controls. A summary of the selected datasets is presented in Table I. All datasets had been generated using the Affymetrix Human Genome U133 Plus 2.0 Array. Among these five studies, one study included two stages of breast cancer (stages I–II), three studies included three stages of breast cancer (stages I–III), and one study included four stages of breast cancer (stages I–IV). As there was only 1 patient in stage IV, the datasets were divided into three subgroups (stage I–III).

R v3.2.2 (https://www.r-project.org/) was used for data pre-processing. The Robust Multichip Average (RMA) algorithm in the oligo-package was used to normalize the raw expression data and generate normalized gene expression intensity (24). Gene annotation, integration and re-normalization of the five datasets were performed using the custom-written Python code (25). Probes with no gene annotation or those that matched multiple gene symbols were removed. Next, the mean expression value of multiple probe IDs that matched an official gene symbol was calculated, and this value was considered to represent the expression intensity of the corresponding gene symbol. The re-normalization method was reported in a previous study (26).

Differential gene expression analysis was performed using R v3.2.2 and the Bioconductor Library (http://www.bioconductor.org/). The empirical Bayes algorithm (function ‘eBayes’) in the limma package was used to detect differentially expressed genes between breast cancer and controls (27). Genes were considered to be upregulated if the logarithmic transformed fold-change log2(FC) was ≥log2(1.5) and the false discovery rate (FDR)-adjusted P-value was ≤0.05. Genes were considered to be downregulated if log2(FC) ≤-log2(1.5) and FDR-adjusted P≤0.05. Differential expression analysis was performed for the whole cohort and the sub-groups (stage I–III). The control samples in the analysis for different stages were the same as the controls in the analysis for the whole cohort.

The Database for Annotation, Visualization and Integrated Discovery Bioinformatics Resources 6.7 was used to perform GO and KEGG pathway enrichment analysis (28). The input parameters were the list of differentially expressed genes. The significance level for enrichment was set at P≤0.05. The 4-set Venn diagram in InteractiVenn (http://www.interactivenn.net/) was used to present the GO terms in the category biological process in which the differentially expressed genes in the unstaged cohort and in different stages of breast cancer were enriched.

Breast cancer-specific target genes were defined based on already available drugs or drugs under development that target these genes. All of these target genes were searched and downloaded from the Thomson Reuters Integrity Database. In total, a list of 344 breast cancer target genes were obtained, which were then mapped to the differentially expressed genes obtained in the present study for the whole cohort and for the sub-groups (stage I–III). The differentially expressed genes identified in the whole cohort and the sub-groups were overlapped with the 344 breast cancer target genes, and the ‘barplot’ function was used to present the results.

A genome-scale integrated analysis of gene networks in breast cancer was performed using the Genome-scale Integrated Analysis of Gene Networks in Tissues (GIANT) web server (http://giant.princeton.edu/) (29). Based on the aforementioned results of the overlapping expression pattern of differentially expressed breast cancer target genes among the whole cohort and the sub-groups, the differentially expressed target genes in the unstaged cohort were used as input parameters to perform the gene network analysis. As the tissue options in the GIANT web server did not include breast tissue, ‘all tissues’ was selected to perform the analysis. The server generated a gene network of target genes and other genes that interacted with the target genes, and biological function enrichment analysis of the genes in the network was performed. The enriched biological processes were then presented using bar charts.

The number of differentially expressed genes in breast cancer for each dataset is presented in Table II. A total of 153 upregulated and 183 downregulated genes were obtained for the whole cohort. In the unstaged cohort and in the stage I–II groups, more downregulated than upregulated genes were identified. However, in the stage-III group, the number of upregulated genes was higher than that of downregulated genes. The number of overlapping up- and downregulated genes among all four groups was 29 and 51, respectively.

A Venn diagram displaying the enrichment results for the GO category biological process for the unstaged cohort and the individual stages is presented in Fig. 1. In the unstaged cohort, stage I–III groups, the differentially expressed genes were enriched in 138, 21, 119 and 136 GO terms, respectively. The top 10 enriched GO terms in the category biological process in each group are presented in Table III. Only one GO term in the category biological process, namely ‘cell migration’, was enriched in all four groups (the P-values for the unstaged cohort, stage I–III groups were 0.018, 0.001, 0.008 and 0.035, respectively). Activated cell migration is known to promote breast cancer progression (30). By contrast, inhibition of breast cancer cell migration contributes to successful treatment (31). Furthermore, 67 GO terms in the category biological process were enriched in the unstaged cohort, stage II–III groups, and 29 were enriched in the unstaged cohort and the stage III group.

The KEGG pathways in which the differentially expressed genes of the four groups were enriched are presented in Table IV. The pathways ‘glutathione metabolism’, ‘peroxisome proliferator activated receptor (PPAR) signaling pathway’, ‘metabolism of xenobiotics by cytochrome P450 (CYP)’, ‘arachidonic acid metabolism’, ‘drug metabolism’ and ‘tight junction’ were significantly enriched in the whole cohort (P-values of <0.001, 0.001, 0.004, 0.018, 0.025 and 0.034, respectively). Fig. 2 presents the gene expression profiles in the sets of enriched pathways in the whole cohort and the three sub-groups based on cancer stage. Overall, all of these pathways were severely affected. The pathways ‘PPAR signaling’, ‘arachidonic acid metabolism’, ‘propanoate metabolism’ and ‘fatty acid metabolism’ had a large number of downregulated genes across all groups. However, the ‘spliceosome’ pathway had more upregulated than downregulated genes. Glutamine has been reported to control cancer cell proliferation by activating signal transducer and activator of transcription 3 independent of glutamine metabolism (32). According to a study on mammary epithelial cell-specific PPARγ knockout mice, PPARγ expression and signaling has an inhibitory role in breast tumor progression (33). According to a previous study by our group, certain downstream genes mainly involved in lipid metabolism and adipocyte differentiation in the PPAR signaling pathway were suppressed following downregulation of PPARγ in breast cancer, which may lead to tumorigenesis (13). It was previously reported that the expression of the CYP1A1, −2E1 and −3A4 was downregulated in tumor tissue, which may alter the biological effects of carcinogens and may represent a potential target for breast cancer chemoprevention (34). Previous studies demonstrated that induction of the expression of CYP2E1 reduces, whereas downregulation of CYP2E1 increases the migratory capacity, thereby promoting breast cancer cell progression (35). Thus, CYP2E1 may be associated with the regulation of breast cancer cell migration. CYP2E1 gene encodes a member of the cytochrome P450 superfamily of enzymes. Cytochrome P450 proteins are monooxygenases that catalyze many reactions involved in drug metabolism and synthesis of cholesterol, steroids and other lipids (36). It was previously demonstrated that differences in the expression of drug and xenobiotic metabolizing enzymes (DXME) markedly affect drug resistance. Substantial differences in DXME expression were identified in breast cancer patients of different ethnicities, which may affect pathways involved in drug metabolism (37).

The differentially expressed breast cancer target genes were screened in the four groups. Subsequently, 8, 6, 9 and 7 breast cancer target genes from the differentially expressed genes in the unstaged cohort and stage I–III groups, respectively, were mapped. The combined set of these target genes contained 14 genes. Fig. 3 presents the log2(FC) of these targets in each group. Insulin-like growth factor (IGF) 1 and TCDD-inducible poly (ADP-ribose) polymerase (TIPARP) were overexpressed in all groups, whereas aryl hydrocarbon receptor nuclear translocator 2 (ARNT2), INSM1, mucin (MUC)1 and PGR were downregulated in all groups, compared with the healthy controls. Overall, all these targets in the groups of different stages exhibited the same expression pattern as in the unstaged cohort. The IGF1 signaling axis has been reported to be crucial for tumorigenesis, and the activation of IGF1 receptor may promote breast cancer development by increasing glycolysis and promoting biomass production (38). Several polymorphisms of IGF1 pathway genes were reported to be associated with the risk of breast cancer (39). TIPARP is a poly(ADP-ribose) polymerase and a transcriptional repressor of the aryl hydrocarbon receptor, the polymorphism of which was previously reported to be associated with ovarian and breast cancer (40,41). The mRNA expression level of ARNT2 was previously reported to be useful in determining the prognosis of breast cancer, and ARNT2 was reported to form functional complexes with hypoxia-inducible factor (HIF), which is a key to factor involved in tumor angiogenesis (42,43). The results of small interfering RNA-mediated knockdown of ARNT2 suggested that ARNT2 may have a pivotal part in the modulation of HIF-1-regulated signaling and metabolism in MCF7 human breast cancer cells (43). The tumor oncoprotein MUC1 is a potential target in breast cancer immunotherapy, and the expression of MUC1 is absent or low in normal breast tissue, while it is high in breast cancer (44). PGR is one of the well-established breast cancer biomarkers, along with HER2/ERBB2 and ER (45).

The present study identified 8 differentially expressed anti-breast cancer target genes in the whole cohort (Fig. 3), which were used to perform a genome-scale integrated analysis. The gene-gene interaction network and the top 15 enriched biological processes are shown in Fig. 4. Among all these targets, RRM2 displayed the highest degree of interaction with other interacting genes. Ribonucleotide reductase M2 (RRM2) is required for pyrimidine metabolism, and it is associated with aggressive tamoxifen-resistant breast tumors, whereas pharmacological inhibition and genetic knockdown of RRM2 sensitizes tumors to tamoxifen (46). In MCF-7 breast cancer cells, overexpression of RRM2 reduced the expression of ERα66 and caused an upregulation of the 36-kDa variant of ER, ERα36, resulting in a reduction in the effectiveness of tamoxifen, which is widely used as an adjuvant therapy for patients with ERα-positive tumors (47). Therefore, RRM2-associated metabolites may potentially be developed as prognostic markers for breast cancer. Furthermore, MUC1 exhibited a high degree of interaction with ERBB2 and ERBB3, PGR exhibited a high degree of interaction with nuclear receptor corepressor 2 and CUE domain containing 2, and cyclin-dependent kinase 5 regulatory subunit associated protein 3 exhibited a high degree of interaction with collapsin response mediator protein 1 (Fig. 4A). It has been demonstrated that the oncogenic MUC1 C-terminal may act on the polycomb repressive complex 1 during epigenetic gene silencing, which is overexpressed in breast and other cancer types (48). The MUC1 oncoprotein was reported to be aberrantly overexpressed and associated with HER2/ERBB2 activation in breast cancer cells (49). The ERBB2/ERBB3 heterocomplex is a vital etiological feature of breast cancer, and it is important to understand its mechanisms of action to improve the design of novel, effective chemotherapeutics (50). As presented in Fig. 4B, the enriched biological processes of these target genes and interacting genes were mostly associated with cell cycle and mitosis. These results indicated that altered expression of these anti-breast cancer genes may severely affect the cell cycle and mitosis. Previous cell cycle-targeting agents have been reviewed, and emerging strategies for targeting mitosis in cancer have been refined and improved (51). Of note, in a recent study on functional mutagenesis screens in mice, human breast cancer susceptibility genes were, at large, not associated with cell cycle/mitosis genes (52). These results suggested that integration of human cancer transcriptomic data is required to identify breast cancer biomarkers with a high prognostic value. A limitation is that the tissue options in the GIANT web server did not include breast tissue, we choose ‘all tissues’ option to perform the analysis. This may have some influence on the result.