Prognostic genes of breast cancer revealed by gene co-expression network analysis.

The aim of the present study was to identify genes that may serve as markers for breast cancer prognosis by constructing a gene co-expression network and mining modules associated with survival. Two gene expression datasets of breast cancer were downloaded from ArrayExpress and genes from these datasets with a coefficient of variation >0.5 were selected and underwent functional enrichment analysis with the Database for Annotation, Visualization and Integration Discovery. Gene co-expression networks were constructed with the WGCNA package in R. Modules were identified from the network via cluster analysis. Cox regression was conducted to analyze survival rates. A total of 2,669 genes were selected, and functional enrichment analysis of them revealed that they were mainly associated with the immune response, cell proliferation, cell differentiation and cell adhesion. Seven modules were identified from the gene co-expression network, one of which was found to be significantly associated with patient survival time. Expression status of 144 genes from this module was used to cluster patient samples into two groups, with a significant difference in survival time revealed between these groups. These genes were involved in the cell cycle and tumor protein p53 signaling pathway. The top 10 hub genes were identified in the module. The findings of the present study could advance the understanding of the molecular pathogenesis of breast cancer.

Introduction

Breast cancer is the most common type of cancer in women, accounting for 25% of all cases (1). Risk factors include lifestyle (including smoking or diet), genetics and medical conditions. A number of treatment methods are now available for breast cancer, including surgery, radiotherapy, chemotherapy, hormone therapy and targeted therapy. However, certain patients have a poor prognosis and the molecular mechanisms underlying this remain unclear. Prognostic factors include disease stage and grade, recurrence of the disease, and the age and health of the patient.

With advances in technology and the accumulation of research results, certain molecular markers associated with breast cancer have been well studied. Tumor protein p53 mutations are poor prognostic factors in breast cancer (2). MYC proto-oncogene and bHLH transcription factor-driven accumulation of 2-hydroxyglutarate are associated with poor breast cancer prognosis (3). Prostaglandin-endoperoxide synthase 2 expression predicts worse breast cancer prognosis (4). Ki-67 has been associated with disease-free survival, but its prognostic value remains to be validated (5). Matrix metalloproteinase-8 gene variation may influence breast cancer prognosis and can have an inhibitory effect on cancer metastasis (6). A gene signature involved in tumor-immune interactions may provide a more accurate prognostic tool (7). Zhang et al (8) performed a meta-analysis and demonstrated that overexpression of C-X-C motif chemokine receptor 4 was significantly associated with lymph node status and distant metastasis, indicating poor overall and disease free survival. SRY-box 4 overexpression is a biomarker for malignant status and poor prognosis in breast cancer patients (9). A number of other novel biomarkers have also been also identified, including chromobox homolog 1 (10), HOX transcript antisense intergenic RNA (9) and anterior gradient 3 (11). Nevertheless, more prognostic genes are required to further improve treatment decisions and thus the quality of life of patients with breast cancer.

Microarray technology has been widely used to identify biomarkers of breast cancer (12,13), allowing for the large-scale screening of molecular markers. In the present study, two gene expression datasets were obtained to reveal prognostic genes (14,15). One dataset was used with the aim of identifying genes associated with the distant metastasis of lymph-node-negative primary breast cancer (14); the other was used to identify genes involved in response and survival following taxane-anthracycline chemotherapy in breast cancer (14). The two datasets were combined to construct a gene co-expression network and analyze survival time to identify novel biomarkers associated with breast cancer prognosis.

Materials and methods

Raw data and pre-treatment

Two gene expression datasets, GSE2034 (14) and GSE25066 (15), were downloaded from ArrayExpress (https://www.ebi.ac.uk/arrayexpress/). Dataset GSE2034 included 286 breast cancer samples and dataset GSE25066 included 508 breast cancer samples. The two gene expression datasets were obtained using Affymetrix GPL96 platform.

Normalization was performed with rma from the affy package (16) in R (R 3.2.0; https://www.r-project.org/) and then log2 conversion was applied. Probes were mapped onto genes according to annotation files. Probes mapping to the same gene were averaged as the expression level for the gene.

Functional enrichment analysis

Gene Ontology (GO) annotation and pathway enrichment analysis were performed with DAVID (Database for Annotation, Visualization and Integration Discovery; http://david.abcc.ncifcrf.gov/) (17).

Gene co-expression network and modules

The gene co-expression network was constructed with the WGCNA package (18) in R. The adjacency coefficient a was calculated as follows:

Where x and x are vectors of expression value for gene i and j; cor represents the Pearson's correlation coefficient of the two vectors; a is the adjacency coefficient and is acquired via exponential transform of S.

WGCNA method takes topological properties into consideration to identify modules from gene co-expression networks. Therefore, this method not only considers the association between the two connected nodes, but also takes associated genes into account. It calculates the weighting coefficient W from a as follows:

W considers the overlap between neighbor genes of genes i and j. Modules were identified via hierarchical clustering of the weighting coefficient matrix, W.

Survival analysis

Cox regression was performed with hub genes from the modules to identify survival-associated genes, and Kaplan-Meier survival was used to compare the survival time of different groups, which were performed with the Survival package in R (https://cran.r-project.org/web/views/Survival.html). P<0.05 was considered to indicate a statistically significant difference. Pearson's correlation was performed by cor function in R (19).

Results

Gene expression data

A total of 13,191 genes were identified in the GSE2034 and GSE25066 datasets, for which box plots are presented in Fig. 1. According to the box plots, the average total mRNA expression level in each sample was consistent, indicating that a good performance of normalization was achieved for both datasets.

Figure 1.

Box plots of normalized gene expression data of two datasets. (A) GSE2034 (286 samples) and (B) GSE25066 (200 samples randomly selected from the total 508 samples). The average total mRNA expression level in each sample was consistent, indicating that a good performance of normalization was achieved. The x-axis represents the gene expression level; the y-axis represents the samples.

A total of 2,669 genes with coefficient of variation (CV) >0.5 were selected. Functional enrichment analysis revealed that they were associated primarily with immune response, cell proliferation, cell differentiation and cell adhesion (Table I).

Table I.

Top 15 significantly over-represented biological pathways.

ID	Description	P-value	Adjusted P-value
GO:0006955	Immune response	2.19×10−63	3.27×10−61
GO:0006952	Defense response	2.15×10−57	2.80×10−55
GO:0006950	Response to stress	1.57×10−56	1.96×10−54
GO:0007166	Cell-surface receptor signaling pathway	1.20×10−55	1.43×10−53
GO:0008283	Cell proliferation	1.06×10−49	1.09×10−47
GO:0002682	Regulation of immune system process	6.12×10−42	4.82×10−40
GO:0016477	Cell migration	7.58×10−40	5.66×10−38
GO:0045321	Leukocyte activation	1.90×10−39	1.32×10−37
GO:0006954	Inflammatory response	3.92×10−38	2.66×10−36
GO:0048584	Positive regulation of response to stimulus	6.10×10−38	4.05×10−36
GO:0042127	Regulation of cell proliferation	1.72×10−37	1.10×10−35
GO:0030154	Cell differentiation	3.01×10−34	1.70×10−32
GO:0048869	Cellular developmental process	2.06×10−33	1.14×10−31
GO:0007155	Cell adhesion	7.77×10−33	4.22×10−31
GO:0022610	Biological adhesion	1.11×10−32	5.90×10−31

Adjusted P-value: Use the multiple comparisons in General Linear Model ANOVA, the adjusted P-value indicates which factor level comparisons within a family of comparisons (hypothesis tests) are significantly different.

Prognostic genes

Two gene co-expression networks were constructed for the two datasets by WGCNA (Fig. 2). Seven modules were identified from the network of GSE2034 via hierarchical clustering of the weighting coefficient matrix, W (Fig. 3). The modules were termed the red, blue, green, black, brown, yellow and turquoise modules.

Figure 2.

Gene co-expression networks for datasets GSE2034 (left) and GSE25066 (right). The x-axis represents the degree of the node, k, while the y-axis represents proportion of genes with degree of k, p (k).

Figure 3.

Seven modules identified from the gene co-expression network. Cluster analysis result is shown above and module identification shown below.

The degree, k, for each gene in the module was calculated and the P-value of Cox regression between each gene and survival was also determined. Next, the correlation between k and -log10 (P) was calculated. The yellow module exhibited significant correlation with survival time in dataset GSE2034 (P=9.3×10−13) (Fig. 4A), which was also observed in dataset GSE25066 (P=9.3×10−6) (Fig. 4B). Besides, survival-associated genes (P<0.05 in Cox regression) were significantly over-represented in the yellow module in both datasets (Fig. 5). Therefore, the yellow module was considered to be significantly associated with breast cancer patient survival, which should be further investigated to understand the association between survival time and critical gene expression.

Figure 4.

Scatter plots of the degree and P-value of Cox regression in datasets (A) GSE2034 and (B) GSE25066. The x-axis indicates the degree of regression, the y-axis indicates the P-value. Each circle represents a gene.

Figure 5.

Survival-associated genes in each module. The x-axis indicates the module, the y-axis indicates the significance of over-representation.

The 144 genes from the yellow module were used in the cluster analysis of samples from dataset GSE2034, which separated the patient samples into two groups based on the expression of these genes (Fig. 6). A significant difference in survival time was observed between the two groups (P=0.008; Fig. 7). Functional enrichment analysis indicated that the 144 genes from the yellow module were involved in cell cycle, oocyte meiosis, the tumor protein p53 signaling pathway and progesterone-mediated oocyte maturation (Table II).

Figure 6.

Cluster analysis using the degree of expression of 144 survival-associated genes for the samples in the GSE2034 dataset.

Figure 7.

Survival curves for the two groups of breast cancer patient samples clustered according to expression of the 144 genes.

Table II.

KEGG pathways enriched in the 144 genes of the yellow module.

ID	Description	P-value	Adjusted P-value
hsa04110	Cell cycle	5.22×10−18	3.13×10−17
hsa04114	Oocyte meiosis	2.17×10−9	6.50×10−9
hsa04115	p53 signaling pathway	2.46×10−5	4.91×10−5
hsa04914	Progesterone-mediated oocyte maturation	9.19×10−5	1.38×10−4

p53, tumor protein p53. Adjusted P-value: Using multiple comparisons in a general linear model analysis of variance, the adjusted P-value indicates which factor level comparisons within a family of comparisons (hypothesis tests) are significantly different.

The top 10 hub genes from the yellow module were selected (Table III) and included cyclin B2 (CCNB2), ubiquitin-conjugating enzyme E2C (UBE2C), protein regulator of cytokinesis 1 (PRC1), cell division cycle 20 (CDC20), abnormal spindle microtubule assembly (ASPM), forkhead box M1 (FOXM1), kinesin family member 4A (KIF4A), nucleolar and spindle associated protein 1 (NUSAP1), pituitary tumor-transforming 1 (PTTG1) and centrosomal protein 55 kDa (CEP55). All of these genes were significantly associated with survival time in the two datasets.

Table III.

Top 10 hub genes in the yellow module.

Dataset	Gene name	Coefficient	P-value	kTotal	kWithin
GSE2034	CCNB2	0.3640	0.0003	14.7998	12.4392
	PRC1	0.3868	0.0005	12.9603	11.3677
	UBE2C	0.4281	0.0006	14.1236	11.2433
	ASPM	0.3442	0.0002	12.9467	10.9328
	CDC20	0.2339	0.0065	14.6847	10.7527
	FOXM1	0.1988	0.0168	13.7352	10.7131
	CEP55	0.3691	0.0004	12.6988	10.6131
	KIF4A	0.2648	0.0217	12.1095	10.3165
	NUSAP1	0.3931	0.0012	11.7988	10.2885
	PTTG1	0.4027	0.0019	12.4981	10.2449
GSE25066	CCNB2	0.323932	0.3239	0.0006	9.4109
	PRC1	0.276034	0.2760	0.0023	6.6109
	UBE2C	0.381925	0.3819	0.0003	6.1036
	ASPM	0.207911	0.2079	0.0031	4.9210
	CDC20	0.329027	0.3290	0.0000	8.5936
	FOXM1	0.170967	0.1710	0.0091	5.9345
	CEP55	0.304415	0.3044	0.0002	6.3694
	KIF4A	0.568168	0.5682	0.0001	3.1945
	NUSAP1	0.270014	0.2700	0.0061	6.7332
	PTTG1	0.791755	0.7918	0.0000	4.0029

Discussion

Two gene expression datasets of breast cancer were obtained and the 2,669 differentially expressed genes with a CV >0.5 were selected. These genes were implicated in the immune response, cell proliferation and cell migration. These functions were closely associated with the development and metastasis of cancer. A breast-cancer-specific gene co-expression network was constructed for dataset GSE2034, from which 7 modules were identified. The yellow module was closely associated with survival time and, as such, the 144 genes from yellow module were investigated further. These genes were primarily involved in the cell cycle and tumor protein p53 signaling pathway. The top 10 hub genes were identified in the yellow module, all of which were associated with poor patient prognosis.

The majority of the 10 critical genes in the yellow module are associated with the cell cycle. CCNB2 is an essential component of the cell-cycle regulatory machinery (20). Elevated CCNB2 expression in invasive breast cancer is associated with unfavorable clinical outcomes (21). UBE2C is required for the degradation of mitotic cyclins and for cell-cycle progression, and is involved in cancer progression. UBE2C is highly expressed in breast microcalcification lesions (22). The prognostic value of UBE2C has been validated in several studies (23–25). microRNA-196a post-transcriptionally upregulates UBE2C and promotes cell proliferation in breast cancer (26). Inhibition of UBE2C reduces proliferation and sensitizes breast cancer cells to radiotherapy and chemotherapy (27), suggesting that it could serve as a potential therapeutic target. CDC20 is a regulatory protein in the cell cycle. Overexpression of CDC20 predicts short-term breast cancer survival (22). ASPM is essential for normal mitotic spindle function and is a marker for vascular invasion, early recurrence and poor prognosis of hepatocellular carcinoma (28). Increased ASPM expression is also associated with enhanced tumor grade and lower survival rates of epithelial ovarian cancer (29). A significant correlation between the expression of the CCNB2 and ASPM proteins is reported (21), which may serve a role in the development of breast cancer.

FOXM1 is a transcriptional activator involved in cell proliferation, which is a downstream target and marker of HER2 overexpression in breast cancer (30). FOXM1 is implicated in the proliferation, migration and invasion of breast cancer cells (31,32) and serves a role in chemotherapy resistance (33,34). KIF4A is an ATP-dependent microtubule-based motor protein that is involved in the intracellular transport of membranous organelles. KIF4A is implicated in doxorubicin-induced apoptosis in breast cancer cells (35). NUSAP1 may be involved in tumorigenesis and in the processes of invasion and progression of breast cancer (36); it influences the DNA damage response by controlling the protein levels of BRCA1 (37). PTTG1 exhibits tumorigenic activity in vivo and is highly expressed in various tumors; it is associated with endocrine therapy resistance in breast cancer (38). PTTG1 may promote tumor malignancy via the epithelial-to-mesenchymal transition and the expansion of the cancer stem cell population (39). CEP55 is also involved in breast cancer progression (40), possibly exerting an oncogenic function via regulation of the phosphoinositide-3 kinase/protein kinase B pathway and midbody fate (41).

PRC1 encodes a protein involved in cytokinesis, specifically the polarization of parallel microtubules, whose expression level changes markedly in the different phases of the cell cycle. PRC1 has been demonstrated to be a substrate of several cyclin-dependent kinases (CDKs); its alternative splicing results in multiple transcript variants (42,43). Although PRC1 serves an important role in the cell cycle, its role in breast cancer remains unclear. The results of the present study indicate that the role of PRC1 in the pathogenesis of breast cancer necessitates further study.

Gene co-expression network analysis revealed several genes of prognostic significance in breast cancer. The majority of these genes have been validated by previous studies; however, the function of certain critical genes identified by gene co-expression network analysis in breast cancer remains unclear, thus providing targets for further studies. These prospective studied may disclose novel biomarkers or provide targets for breast cancer therapies.