Comprehensive integrated analysis of gene expression datasets identifies key anti-cancer targets in different stages of breast cancer.

Breast cancer is one of the primary threats to women's health worldwide. However, the molecular mechanisms underlying the development of breast cancer remain to be fully elucidated. The present study aimed to investigate specific target gene expression profiles in breast cancer tissues in general and in different breast cancer stages, as well as to explore their functions in tumor development. For integrated analysis, a total of 5 gene expression profiling datasets for 3 different stages of breast cancer (stages I-III) were downloaded from the Gene Expression Omnibus of the National Center for Biotechnology Information. Pre-processing of these datasets was performed using the Robust Multi-array Average algorithm and global renormalization was performed for all studies. Differentially expressed genes between breast cancer patients and controls were estimated using the empirical Bayes algorithm. The Database for Annotation, Visualization and Integrated Discovery web server was used for analyzing the enrichment of the differentially expressed genes in Gene Ontology terms of the category biological process and in Kyoto Encyclopedia of Genes and Genomes pathways. Furthermore, breast cancer target genes were downloaded from the Thomson Reuters Integrity Database. We merged these target genes with the genes in breast cancer datasets. Analysis of anti-breast cancer gene networks was performed using the Genome-scale Integrated Analysis of Gene Networks in Tissues web server. The results demonstrated that the normal functions of the cell cycle, cell migration and cell adhesion were altered in all stages of breast cancer. Furthermore, 12 anti-breast cancer genes were identified to be dysregulated in at least one of the three stages. Among all of these genes, ribonucleotide reductase regulatory subunit M2 (RRM2) exhibited the highest degree of interaction with other interacting genes. Analysis of the network interactions revealed that the transcription factor of RRM2 is crucial for cancer development. Other genes, including mucin 1, progesterone receptor and cyclin-dependent kinase 5 regulatory subunit associated protein 3, also exhibited a high degree of interaction with the associated genes. In conclusion, several key anti-breast cancer genes identified in the present study are mainly associated with the regulation of the cell cycle, cell migration, cell adhesion and other cancer-associated cell functions, particularly RRM2.

Introduction

Breast cancer is one of the most common cancer types among women worldwide. According to global cancer statistics, an estimated >1.6 million patients were newly diagnosed and 500,000 breast cancer-associated mortalities occurred in 2012 worldwide (1). The Global Burden of Disease estimated that there were more than 1.7 million new cases and more than 545 thousand deaths in 2016 (2). In China, breast cancer alone is estimated to account for 15% of all newly diagnosed cancers in women, and its incidence has increased in the past decades (3). At present, chemotherapy is an important means of systemic therapy for breast cancer, in addition to surgical treatment. However, breast cancer may still be associated with poor prognosis, short survival time and rapid recurrence after therapy (4). The treatment and prognosis of breast cancer are affected by the expression levels of certain genes and proteins. For instance, triple-negative breast cancer [estrogen receptor (ER)-, progesterone receptor (PGR)- and human epidermal growth factor receptor 2 (HER2)-negative] is associated with poor prognosis and no targeted systemic therapy is currently available (5). However, patients with triple-positive breast cancer have a better prognosis and longer overall survival compared with triple-negative breast cancer patients (6).

The American Joint Committee on Cancer (AJCC) breast cancer staging system provides important information for the treatment and prognosis of this type of cancer (7). According to the 8th edition of the AJCC staging system, breast cancer may be divided into four main stages (stages I–IV) based on various factors, including the size of the tumor, the status of the lymph nodes and metastasis (7). Breast cancer staging is crucial for determining the extent of disease progression, as well as for containing and eliminating the cancer. The treatment of breast cancer depends partly on the stage of the disease, particularly in the case of targeted therapy. During breast cancer progression, diverse genetic signatures have been identified to drive processes of genome, transcriptome and epigenome remodeling (8–13). Therefore, it is necessary to select the most effective treatment options for breast cancer patients at different stages.

To date, numerous genes and pathways have been identified to be associated with breast cancer, and this information may be useful for studies into the pathological mechanisms and clinical treatment of breast cancer. According to the evidence provided by a functional study, the Wnt/β-catenin signaling pathway controls cell fate in developmental processes and tumorigenesis, with β-catenin identified as a transforming factor (14). Based on a genome-wide assessment of allelic imbalances, an ATR/ATM-regulated DNA damage response network was identified to be activated in early human tumorigenesis, which may delay or prevent tumor progression (15). The results of a medical genomics study indicated that paired-box gene 2 may mediate endometrial carcinogenesis induced by tamoxifen, which has been widely used in the treatment of hormone-responsive breast cancer at all stages (16).

Thomson Reuters Integrity™ is a knowledge solution integrating biology, chemistry and pharmacology data (https://thomsonreutersintegrity.com). It contains exhaustive information on therapeutic drugs and gene targets for numerous complex human diseases. To the best of our knowledge, no previous study has compared the expression profiles of breast cancer target genes at different tumor stages. Gene expression profiling datasets at different stages of breast cancer are available from public databases, including the Gene Expression Omnibus of the National Center for Biotechnology Information (NCBI-GEO; http://www.ncbi.nlm.nih.gov/geo) (17,18). The present study aimed to perform a comprehensive integrated analysis of gene expression datasets to identify key targets for breast cancer treatment and explore similarities and differences in the abnormalities of molecular signaling pathways/biological functions at different stages of breast cancer.

Materials and methods

Microarray data collection and pre-processing

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).

Table I.

Summary of Gene Expression Omnibus breast cancer datasets used in the present study.

Dataset ID	Author (year)	Samples (n)	Breast cancer stages	(Refs.)
GSE10810	Pedraza (2010)	58	I, II, III	(19)
GSE16391	Desmedt (2009)	48[a]	I, II	(20)
GSE29431	Lopez (2012)	66	I, II, III	(21)
GSE42568	Clarke (2013)	121	I, II, III	(22)
GSE61304	Aswad (2015)	62	I, II, III, IV	(23)

The dataset contained 55 samples, but only 48 were used, as the others lacked case/control information. In the other datasets, all samples were used.

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

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.

Enrichment analysis for differentially expressed genes in Gene Ontology (GO) terms in the category biological process and in KEGG pathways

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 target gene analysis

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.

Gene interaction network analysis

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.

Results and Discussion

Overview of differentially expressed genes

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.

Table II.

Number of differentially expressed genes in breast cancer.

Group	Cases[a]	Mapped genes	Upregulated	Downregulated
Entire cohort	257	20307	153	183
Stage I	22	20307	53	275
Stage II	98	20307	167	309
Stage III	113	20307	202	165

Values are expressed as n.

With 98 controls in each group.

GO and KEGG enrichment results

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.

Figure 1.

Venn Diagram of enriched Gene Ontology terms in the category Biological Process in breast cancer. The four groups (unstaged cohort, stage I–III) are represented by red, blue, cyan and orange color, respectively. The unstaged cohort contained 1 stage IV sample and 23 samples without stage information.

Table III.

Top 10 enriched GO terms in the category biological process by the differentially expressed genes from the gene expression datasets for breast cancer.

Group/GO term	P-value
Entire cohort
Response to wounding	<0.001
Epithelial cell differentiation	<0.001
Response to endogenous stimulus	<0.001
Response to nutrient levels	<0.001
Epithelium development	<0.001
Regulation of hormone levels	<0.001
Defense response	<0.001
Response to drug	<0.001
Response to extracellular stimulus	<0.001
Response to steroid hormone stimulus	<0.001
Stage I
Cell migration	0.001
Vasculature development	0.002
Localization of cell	0.004
Cell motility	0.004
Endothelial cell migration	0.005
Angiogenesis	0.009
Odontogenesis	0.010
Leukocyte migration	0.012
Blood vessel development	0.016
Blood vessel morphogenesis	0.019
Stage II
Gland development	<0.001
Response to extracellular stimulus	<0.001
Cellular di-, tri-valent inorganic cation homeostasis	<0.001
Response to wounding	<0.001
Di-, tri-valent inorganic cation homeostasis	<0.001
Response to nutrient levels	<0.001
Response to nutrient	<0.001
Cellular cation homeostasis	<0.001
Cell-cell signaling	<0.001
Regulation of hormone levels	<0.001
Stage III
Response to endogenous stimulus	<0.001
Response to hormone stimulus	<0.001
Response to steroid hormone stimulus	<0.001
Response to organic substance	<0.001
Response to nutrient levels	<0.001
Response to wounding	<0.001
Oxidation reduction	<0.001
Epithelial cell differentiation	<0.001
Defense response	<0.001
Response to oxygen levels	<0.001

GO, gene ontology.

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).

Table IV.

Enriched KEGG pathways by the differentially expressed genes from the gene expression datasets for breast cancer.

Group/KEGG pathway	P-value
Entire cohort
Glutathione metabolism	<0.001
PPAR signaling pathway	0.001
Metabolism of xenobiotics by cytochrome P450	0.004
Arachidonic acid metabolism	0.018
Drug metabolism	0.025
Tight junction	0.034
Stage I
Small cell lung cancer	0.008
Focal adhesion	0.009
ECM-receptor interaction	0.035
Spliceosome	0.038
Cytosolic DNA-sensing pathway	0.047
Stage II
Tight junction	0.034
Focal adhesion	0.041
ECM-receptor interaction	0.043
Vascular smooth muscle contraction	0.043
Metabolism of xenobiotics by cytochrome P450	0.048
Stage III
Propanoate metabolism	0.001
Glutathione metabolism	0.001
PPAR signaling pathway	0.001
Metabolism of xenobiotics by cytochrome P450	0.002
Fatty acid metabolism	0.011
Arachidonic acid metabolism	0.033
Glycolysis/gluconeogenesis	0.041
Drug metabolism	0.046

KEGG, Kyoto Encyclopedia of Genes and Genomes; ECM, extracellular matrix; PPAR, peroxisome proliferator activated receptor.

Figure 2.

Gene expression profiles of enriched KEGG pathways in breast cancer. The 15 enriched KEGG pathways are represented by different colors. The red bars represent the upregulated genes and the blue bars represent the downregulated genes. KEGG, Kyoto Encyclopedia of Genes and Genomes; ECM, extracellular matrix; PPAR, peroxisome proliferator activated receptor. The unstaged cohort contained 1 stage IV sample and 23 samples without stage information.

Mapping of anti-breast cancer target genes

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).

Figure 3.

LogFC bar graph of mapped breast cancer-associated genes. (A) LogFC in the unstaged cohort. (B) LogFC in different stage groups. The horizontal dashed lines represent the logFC cut-off for the up- and downregulated genes. *False discovery rate-adjusted P<0.05. FC, fold change.

Gene network of breast cancer targets

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.

Figure 4.

Genome-scale integrated analysis of gene networks in breast cancer. (A) The gene interaction network of breast cancer-associated target genes (query genes) and associated genes (other genes). (B) Top 15 Gene Ontology terms in the category biological process among the breast cancer-associated genes in the network. FDR, false discovery rate.