Prognostic Genes of Breast Cancer Identified by Gene Co-expression Network Analysis

1Breast cancer is the most frequently diagnosed malignancy and the second leading cause of cancer death in females worldwide, accounting for 30% of cancer diagnoses and 14% of cancer death. In 2017, it was estimated that nearly 252,710 new cases were diagnosed in the United States, with ~40,610 deaths. Therapeutic strategies of breast cancer have been markedly improved. A number of treatments such as surgery, chemotherapy, radiotherapy, hormone therapy, and targeted therapy are available for breast cancer. However, the patients with distant metastases were usually diagnosed with a late stage and nearly incurable. Moreover, 30% patients diagnosed with early stage were easy to recur in distant organs even after surgery of removing the primary tumor. The classification of breast cancer affects treatment decision and prognosis: hormone-based therapy for ER+ patients; targeted therapy for HER2+ patients; and poorly differentiated cancer often has the worse prognosis.

2Inheritance plays an important role in the development of breast cancer. BRCA1 and BRCA2 are 2 biomarkers which are currently used clinically to assess the familial breast cancer risk. BRCA-associated breast cancer has relatively distinct pathologic characteristics. Up to 20% women with triple-negative breast cancer present BRCA mutations, while BRCA mutations occur less common in general population. HER2 expression was found to be upregulated in over 30% patients with breast cancer. Previous data suggested that high HER2 levels not only indicated prognostic value, but also affected treatment decisions. Lapatinib and trastuzumab presented dramatically therapeutic effects in patients with HER2-positive breast cancer. Expression levels of hormone receptors (ER/PR) predicted the efficacy of endocrine therapies, and their upregulation was often associated with a favorable prognosis. Ki-67 was reported to be associated with disease-free survival. High CXCR4 levels were associated with lymph node metastasis and distant metastasis. Despite the substantial improvements in the treatment of breast cancer, to date, the ability to treat the advanced ones is still limited due to the lack of precise molecular targets for breast cancer. Therefore, it is important to explore the molecule mechanisms involved in the occurrence and development of breast cancer. More novel candidate genes are needed to improve the early diagnosis and treatment decisions.

3Co-expression analysis is a powerful technique to construct free-scale gene co-expression networks. The weighted gene co-expression network analysis (WGCNA) was widely used to analyze large-scale data sets and to find modules of highly correlated genes. WGCNA was successfully used to explore the associations between gene sets and clinical features, and to identify candidate biomarkers. Thus, we described the correlation patterns among genes through a systematic biology method based on WGCNA and identified novel biomarkers associated with breast cancer prognosis.

4A workflow of this study was indicated in Figure 1. The gene expression profiles of GSE1561 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE1561) submitted by Richard Iggo et al. was downloaded from the Gene Expression Omnibus (GEO) database. The GSE1561 was an expression profiling based on GPL96 platform (Affymetrix Human Genome U133A Array) and contained 49 samples. Most patients had 2 trucut biopsies taken, and both biopsies were analyzed from 2 tumors to test the reproducibility of the technique. Repeat amplifications and duplicate biopsies clustered together suggested that biological variation was greater than technical variation in this data set. The results of immunohistochemistry (IHC) also suggested the high quality of this data set. Robust Multi-array Average (RMA) algorithm in affy package within Bioconductor (http://www.bioconductor.org) in R was used to preprocess the gene expression profile data. After background correction, quantile normalization and probe summarization, the data set with 12,413 genes was further processed, and the top 50% most variant genes by analysis of variance (6,206 genes) were selected for WGCNA analysis.

5After validation, the expression data profile of these 6,206 genes were constructed to a gene co-expression network using WGCNA package in R (Supplementary Data Sheet 1). The analysis was performed as described previously.

6The adjacency matrix aij which calculated the connection strength between each pair of nodes was calculated as follows: Where Xi and Xj were vectors of expression value for gene i and j, sij represented the Pearson's correlation coefficient of gene i and gene j, aij encoded the network connection strength between gene i and gene j. In the presented study, the power of β = 9 (scale free R2 = 0.95) was selected as the soft-thresholding parameter to ensure a scale-free network. In the co-expression network, genes with high absolute correlations were clustered into the same module. WGCNA method not only considers the association between the 2 connected genes, but also takes associated genes into account. Modules were also identified via hierarchical clustering of the weighting coefficient matrix. To further identify functional modules in the co-expression network with these 6,206 genes, the topological overlap measure (TOM) representing the overlap in shared neighbors, was calculated using the adjacency matrix. Where A is the weighted adjacency matrix given by and is the soft thresholding power. According to the TOM-based dissimilarity measure with a minimum size (gene group) of 30 for the gene dendrogram, average linkage hierarchical clustering was conducted, and genes with similar expression profiles were classified into the same gene modules using the DynamicTreeCut algorithm.

7Two approaches were used to identify modules associated with clinical information of breast cancer. First, module eigengenes (MEs) were defined as the first principal component of each gene module and the expression of MEs was considered as a representative of all genes in a given module. The correlation between MEs and clinical trait was calculated to identify the clinical significant module. In addition, the gene significance (GS) was defined as mediated p-value of each gene (GS = lgP) in the linear regression between gene expression and the clinical traits. Then, the module significance (MS) were defined as the average GS of all the genes involved in the module.MS was measured to incorporate clinical information into the co-expression network. Module significance (MS) was defined as the average absolute gene significance measured for all genes in a given module.

8DAVID (http://david.abcc.ncifcrf.gov/) is a database for annotation, visualization and integrated discovery. Gene Ontology (GO) and KEGG pathway analysis of differentially expressed mRNAs were carried out using DAVID (version 6.8) online tools: functional annotation. The ontology contains three categories: biological process (BP), molecular function (MF), and cellular component (CC). Enriched GO terms and KEGG pathways were identified according to the cut-off criterion of adjusted P < 0.001.

9The connectivity of genes was measured by absolute value of the Pearson's correlation. Genes with high within-module connectivity were considered as hub genes of the modules (cor.geneModuleMembership > 0.8). Hub genes inside a given module tended to have a strong correlation with certain clinical trait, which was measured by absolute value of the Pearson's correlation (cor.geneTraitSignificance > 0.2). To validate the hub genes, the clinical information and RNA sequencing data of breast cancer were obtained from the Cancer Genome Atlas Project database (TCGA, https://cancergenome.nih.gov/). The mRNA sequencing data was normalized using edgeR package in R language. The Human Protein Atlas (http://www.proteinatlas.org) was also used to validate the immunohistochemistry of candidate hub genes. The direct link to these images in the human protein atlas are as follows: http://www.proteinatlas.org/ENSG00000112029-FBXO5/tissue/breast#img (FBXO5 in normal tissue); http://www.proteinatlas.org/ENSG00000112029-FBXO5/pathology/tissue/breast$+$cancer#img (FBXO5 in tumor tissue); http://www.proteinatlas.org/ENSG00000157456-CCNB2/tissue/breast#img (CCNB2 in normal tissue); http://www.proteinatlas.org/ENSG00000157456-CCNB2/pathology/tissue/breast$+$cancer#img (CCNB2 in tumor tissue); http://www.proteinatlas.org/ENSG00000090889-KIF4A/tissue/breast#img (CCNB2 in normal tissue); http://www.proteinatlas.org/ENSG00000090889-KIF4A/pathology/tissue/breast$+$cancer#img (CCNB2 in tumor tissue); http://www.proteinatlas.org/ENSG00000065328-MCM10/tissue/breast#img (MCM10 in normal tissue); http://www.proteinatlas.org/ENSG00000065328-MCM10/pathology/tissue/breast$+$cancer#img (MCM10 in tumor tissue); http://www.proteinatlas.org/ENSG00000088325-TPX2/tissue/breast#img (TPX2 in normal tissue); http://www.proteinatlas.org/ENSG00000088325-TPX2/pathology/tissue/breast$+$cancer#img (TPX2 in tumor tissue). Survival analysis of hub genes were performed using Kaplan Meier-plotter (www.kmplot.com).

10The samples of GSE1561 were clustered using average linkage method and Pearson's correlation method (Figure 2). The co-expression analysis was carried out to construct the co-expression network. In this study, the power of β = 9 (scale free R2 = 0.95) was selected as the soft-thresholding parameter to ensure a scale-free network (Figure 3). A total of 18 modules were identified via the average linkage hierarchical clustering. Blue module was found to have the highest association with tumor grade (Figure 4), and this module was selected as the clinical significant module for further analysis.

11The genes in the clinical significant module were categorized into 3 functional groups (BP, CC, and MF). Clinical significant module genes in the BP group were mainly enriched in cell division, DNA replication, sister chromatid cohesion, mitotic nuclear division, and DNA replication initiation; The genes in the MF group were mainly enriched in protein binding, poly(A) RNA binding, RNA binding, and ATP binding; the genes in the CC group were significantly enriched in nucleoplasm, nucleus, nucleolus, cytosol, and cytoplasm (Figure 5). According to Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, our results demonstrated that these genes were mainly involved in cell cycle, DNA replication, spliceosome, ribosome biogenesis in eukaryotes and RNA transport. These results indicated that the clinical significant module genes were mainly involved in mitotic cell cycle process.

12Based the cut-off criteria (|MM| > 0.8 and |GS| > 0.2), 42 genes with high connectivity in the clinical significant module were identified as hub genes. Among them, CCNB2, FBXO5, KIF4A, MCM10, and TPX2 were negatively associated with the overall survival and relapse free survival (Figures 6, 7). Moreover, based on the TCGA data, the expression levels of these 5 genes were significantly higher in tumor tissues, especially in the triple negative breast cancers. The expression of CCNB2, KIF4A, and TPX2 were upregulated in the advanced tumor stages. ROC curve indicated that CCNB2, FBXO5, KIF4A, MCM10, and TPX2 exhibited excellent diagnostic efficiency for normal and tumor tissues (Figures 8, 9). In addition, the protein levels of these 5 genes were significantly higher in tumor tissues compared with normal tissues based on the Human Protein Atlas database (Figure 10). Since these 5 genes were all hub genes in the clinical significant module, they might have a tendency to co-express. Our results of correlation analysis demonstrated a strong correlation of mRNA expression levels between KIF4A and TPX2 (Supplementary Data Sheet 2).

13Breast cancer seriously endangers female health, and it is easy to recur even after combined therapy. Although the treatment of breast cancer was improved during the last decades, the ability to treat the advanced ones is still limited due to the lack of precise molecular targets for breast cancer. Therefore, it is important to explore the molecule mechanisms involved in the occurrence and development of breast cancer. Better biomarkers for cancer specific prognosis and progression are highly demanded. In the presented study, we used gene expression datasets from GEO database to screen potential biomarkers related to the progression and prognosis of breast cancer. We also obtained the clinical information and RNA sequencing data of breast cancer from TCGA database for validation.

14WGCNA was performed to explore gene co-expression modules associated with progression of breast cancer. A total of 6,206 most variant genes were used to construct co-expression network and 18 modules were identified. Blue module was found to have the highest association with tumor grades and 42 genes with high connectivity were screened out from the module. Among them, CCNB2, FBXO5, KIF4A, MCM10, and TPX2 were negatively associated with the overall survival (Figure 6).

15CCNB2, also known as cyclin B2, is a member of cyclin family. CCNB2 was reported to regulate cell cycle by activating CDC2 kinase in eukaryotes, and inhibition of CCNB2 induced cell cycle arrest. CCNB2 was overexpressed in multiple tumors, including bladder cancer, uterine corpus endometrial carcinoma, prostate cancer, and gastric cancer. In addition, compared with normal controls, the levels of serum circulating CCNB2 are higher in digestive tract cancer and lung cancer patients, and they are found to be significantly associated with tumor stage and metastasis status. In invasive breast carcinoma, cytoplasmic CCNB2 protein levels were significantly correlated with a poor disease specific survival. CCNB2 expression level was reported to be an independent prognostic factor for the disease specific survival of breast cancer. Our results indicated that CCNB2 was upregulated in breast cancer tissues compared to normal tissues, and that its expression was significantly associated with molecular subtypes of breast cancer and tumor stages (Figure 8). The underlying mechanisms of CCNB2 on tumor progression need to be further clarified.

16F-Box Protein 5 (FBXO5) is a key cell cycle regulatory gene which regulates the progression to S phase and mitosis by inhibiting the anaphase promoting complex (APC). FBXO5 is overexpressed in various solid tumors. In the G0 and early G1 phases, the expression of FBXO5 is low, while in the S phase it is upregulated. In ovarian clear cell carcinoma, FBXO5 accumulation was related to mitotic errors with centrosome overduplication and abnormal spindle formation. These findings demonstrated that it might be involved in human cell cycle disorders and genomic stability to promote tumor growth. In breast carcinoma tissues, FBXO5 induced proliferation through the PI3K/Akt pathway. Overexpression of FBXO5 was reported to correlate with poor prognosis. In addition, PI3K inhibitor reduced FBXO5 expression.

17The protein encoded by Kinesin family member 4A (KIF4A) was reported to be involved in the intracellular transport of membranous organelles and chromosome integrity during mitosis. In patients with colorectal cancer, KIF4A was upregulated, and downregulation of KIF4A reduced cell proliferation in colorectal cancer cells. In hepatocellular carcinoma (HCC) patients, KIF4A overexpression was associated with poorer overall and disease-free survival. In HCC cells, higher levels of KIF4A dramatically increased cellular clonogenic abilities and proliferation, while KIF4A depletion caused a significant augmentation of apoptosis. In breast cancer, high KIF4A levels were associated with poor relapse-free survival of ER-positive patients. In tamoxifen-resistant and sensitive breast cancer cells, KIF4A knockdown significantly impeded cellular proliferation and induced apoptosis.

18Mini-chromosome maintenance complex component 10 (MCM10) is one of the highly conserved mini-chromosome maintenance proteins. MCM10 is bound to chromatin through the interaction with MCM2-7, and plays crucial roles both in initiation and elongation during eukaryotic genome replication. For urothelial carcinoma, high MCM10 levels were significantly correlated with advanced tumors stages, vascular invasion, and nodal status. MCM10 overexpression also predicted poor disease-specific survival and inferior metastasis-free survival. In our analysis of GSE1561, MCM10 was one of the hub genes in the blue module which was significantly associated with tumor grade (Figure 3). In the validation dataset of TCGA, our results indicated that MCM10 was significantly upregulated in breast tumor tissues, and even higher in the triple negative breast cancer (Figures 8, 9).

19Targeting protein for Xenopus kinesin-like protein 2 (TPX2) plays a critical role in chromosome segregation machinery during mitosis. It was reported to be overexpressed in multiple tumors: lung cancer, kidney renal clear cell carcinoma, hepatocellular Carcinoma, prostate cancer, and breast cancer. TPX2 activates PI3K/Akt pathway and upregulates matrix metalloproteinases (MMP) family members in colon cancer. Previous studies showed that TPX2 expression promoted proliferation, migration, and invasion of liver cancer and breast cancer cells via upregulating expressions of MMP2 and MMP9. In patients with HCC, overexpression TPX2 was correlated with worse prognosis. In addition, knockdown TPX2 in HCC cells strongly reduced cellular proliferation, induced apoptosis and inhibited EMT.

20Co-expression analysis is a powerful technique for multigene analysis of large-scale data sets. In cancer research, co-expression analyses revealed the mRNA and microRNA expression network in multiple cancers. In the present study, we used WGCNA to construct a gene co-expression network, to measure the relationships between genes and modules, and to explore the relationships between modules and clinical traits. We also screened out a clinical significant module which was associated with the progression of breast cancer. KEGG pathway analysis demonstrated that this module was mostly involved in cell cycle. In addition, 5 hub genes, CCNB2, FBXO5, KIF4A, MCM10, and TPX2 were identified and validated to be associated with the progression and worse prognosis of breast cancer. Our results provided valuable indication for basic and clinical research on breast cancer. The underlying concept of gene co-expression analysis is guilt-by-association. The groups of genes known as co-expression modules were found to maintain a consistent expression relationship independent of phenotype, and might share a common biological role. Similar to the limitations of most other data mining methods, our results of WGCNA can be biased or invalid when dealing with technical artifacts or tissue contaminations. To increase the credibility of WGCNA results, TCGA RNA-seq data and IHC data from the Human Protein Atlas database were used for validation. While due to the limitation of the database, the related IHC of each sample can't be found, tumor and normal samples were from different patients.

21JT, YG, and GW reviewed relevant literature and drafted the manuscript. DK, KW, DZ, and QC conducted all statistical analyses. All authors read and approved the final manuscript.

22The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.