Comprehensive analysis of biomarkers for prostate cancer based on weighted gene co-expression network analysis.

BACKGROUND: Prostate cancer (PCa) is one of the leading causes of cancer-related death. In the present research, we adopted a comprehensive bioinformatics method to identify some biomarkers associated with the tumor progression and prognosis of PCa.
METHODS: Differentially expressed genes (DEGs) analysis and weighted gene co-expression network analysis (WGCNA) were applied for exploring gene modules correlative with tumor progression and prognosis of PCa. Clinically Significant Modules were distinguished, and Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were used to Annotation, Visualization and Integrated Discovery (DAVID). Protein-protein interaction (PPI) networks were used in selecting potential hub genes. RNA-Seq data and clinical materials of prostate cancer from The Cancer Genome Atlas (TCGA) database were used for the identification and validation of hub genes. The significance of these genes was confirmed via survival analysis and immunohistochemistry.
RESULTS: 2688 DEGs were filtered. Weighted gene co-expression network was constructed, and DEGs were divided into 6 modules. Two modules were selected as hub modules which were highly associated with the tumor grades. Functional enrichment analysis was performed on genes in hub modules. Thirteen hub genes in these hub modules were identified through PPT networks. Based on TCGA data, 4 of them (CCNB1, TTK, CNN1, and ACTG2) were correlated with prognosis. The protein levels of CCNB1, TTK, and ACTG2 had a degree of differences between tumor tissues and normal tissues.
CONCLUSION: Four hub genes were identified as candidate biomarkers and potential therapeutic targets for further studies of exploring molecular mechanisms and individual therapy on PCa.

Introduction

Prostate cancer (PCa) is the foremost regularly non-cutaneous malignancy among males and the second primary cause of cancer-related deaths in the Western world, the occurrence rate of which keeps increasing.[ Although plenty of PCas are characterized as slow-growing and indolent, some PCa patients are at high risk of developing local recurrence and/or distant metastases and die of their disease. The universal diagnosis approaches of PCa contain digital examination of the rectum, prostate-specific antigen (PSA) test, and subsequential biopsies for histopathological staging.[ It is of great benefit but also huge challenges to make the stratification of high-risk PCa patients after diagnosing. At present, management modes of PCa are principally guided by risk degrees on the strength of Gleason grade, clinical stage, as well as serum PSA level. A better understanding of the potential mechanisms and molecular pathways contributing to tumorigenesis and tumor progression is quite essential to enhance our current diagnostic, prognostic, and therapeutic abilities.

Benefited from the development of high-throughput technology, such as whole genome[ and whole exome sequencing,[ researchers could explore and discover the full molecular landscapes of tumors at a variety of levels including DNA structural changes in the PCa genome. In recent years, gene expression profiles, noncoding RNA profiles and DNA methylation profiles are promising not only for serving as prognosis biomarkers but also therapy targets. For example, cell cycle genes are proved as potential markers that could predict the risk of clinicopathological outcomes such as biochemical recurrence (BCR) rate after radical prostatectomy.[ More recently, new research figured out that up-regulated KIF20A expression level was significantly associated with higher BCR, serum PSA value, Gleason score, pathological T stage, and lymph node metastasis in PCa.[ Besides, the down-regulation of miR-24 was reported to promote cell cycle, proliferation, migration, and clonogenicity potential of PCa cells through regulating the expression of p16 and p27, which supported that it may act as a tumor suppressor in PCa.[ Particularly, in human prostate adenocarcinomas, 128 CpG sites were found to have high accuracy and sensitivity in differentiating samples, which could help us identify lymphatic metastases at DNA methylation level and also may inspire novel therapeutic methods.[ Thus, there is a significant clinical need for finding and identifying united novel biomarkers related to clinicopathological parameters, guiding the surveillance of tumorous progression, assessment of prognosis and improvement of personalized treatment programs, especially at gene expression status.

Nevertheless, many researches to date just focused on the screening of differentially expressed genes, and the high degree of interconnection among genes has been largely neglected, where genes with semblable expression patterns might be functionally correlated. Weighted gene co-expression network analysis (WGCNA), a capable approach widely applied for constructing free-scale gene co-expression networks and exploring the connections between different gene sets or between gene sets and interesting phenotypic traits.[ WGCNA has been used successfully in various biological processes, which vastly favorable for the identification of candidate biomarkers or potential treatment targets.[ Moreover, it is of great success to identify biomarkers for castration-resistant prostate cancer by using WGCNA.[ Thus, in the current research, we attempt to construct a free-scale gene co-expression network of correlation patterns among genes through a comprehensive bioinformatics method in the light of WGCNA and identified network-centric genes relative to clinical traits of PCa.

Materials and methods

Microarray data

Figure S1 showed the workflow of our research. We downloaded the gene expression profiles of GSE69223[ (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE69223) from the Gene Expression Omnibus (GEO) database. Data set GSE69223, performed on Affymetrix Human Genome U133 Plus 2.0 Array, contained malignant and matched normal adjacent tissue samples of 15 prostate cancer patients with prostatectomy. It included microarray data of these patients and detailed clinical variables, such as patients age, tumor stage, and grade, which could be used for constructing co-expression networks, identifying clinical features related modules, and screening hub genes. Ethical approval is not necessary because the data of this study is derived from gene database.

Data preprocessing and screening of DEGs

The gene expression profile data were preprocessed through Robust Multi-array Average (RMA) algorithm in affy package within Bioconductor (http://www.bioconductor.org) in R, and quantile normalization was used to normalize their processed signals. Affymetrix annotation files were applied for annotating the probes. If different probes linked to the same gene, the median value was chosen as the gene expression value. DEGs between prostate cancer and matched normal adjacent tissue were analyzed by the “limma” R package in R. |Log2-fold change (LFC)| >0.5 and Pfalsediscoveryrate(FDR) < 0.05 and were considered as the thresholds for DEGs screening.

Construction of co-expression network

Firstly, DEGs were used to construct scale-free co-expression networks. Secondly, the correlation matrix between 2 genes was performed using Pearson correlation matrices and the major connecting rod. Next, the correlation matrix was transformed into a weighted adjacency matrix through a power function axy=|cxy|β (amn=adjacency between gene x and gene y; cmn=Pearson correlation between gene x and gene y). β was the most appropriate soft-thresholding parameter which could raise the co-expression similarity and achieve consistent scale-free topology in multiple datasets. Here, we chose the power of β = 10 (scale-free R2 = 0.85) to ensure a scale-free network (Fig. S2). After that, with the chosen power value, we performed automatic network construction and module detection with the following major parameters: maxBlockSize of 4000, minModuleSize of 100, deepSplit of 2, and mergeCutHeight of 0.5. This procedure comprised the calculation of network adjacencies and topological overlap dissimilarities, followed by scaling of topological overlap matrices (TOM) and calculation of consensus topological overlap. And then, we established an average linkage hierarchical clustering dendrogram of gene expression data for the dataset to cluster genes that have similar expression profiles into 1 gene module and performed adaptive branch cutting to identify shared functional modules. The Above steps were implemented by using the “WGCNA” package in R.

Identification of clinically significant modules

Module eigengenes (MEs) were characterized as the first principal component of each co-expression module and the expression of MEs was regarded as a representative of entire genes in a given module. Additionally, we recognize the clinically significant module (hub module) by calculating the correlation between MEs and clinical data.

Gene ontology and KEGG enrichment analysis

For further understanding the function of DEGs in the hub module, all genes in the hub module were subjected to GO and KEGG pathway analysis through the Database for DAVID database (version 6.8) (http://david.abcc.ncifcrf.gov/). The ontology includes 3 categories: biological process (BP), cellular component (CC), and molecular function (MF). P value < .05 was regarded as the cut-off criterion for identifying the enriched GO terms and KEGG pathways.

Hub gene analysis and construction of PPI network

In this study, after the hub module was confirmed, we selected the candidate hub genes by module connectivity (cor.geneModuleMembership (MM) > 0.6 or >0.8) and clinical feature relationship (cor.geneTraitSignificance (GS) > 0.2), both measured by absolute value of the Pearson correlation. After that, we uploaded these genes to the STRING database (https://string-db.org/) for constructing protein–protein interaction (PPI) networks and screening out hub nodes. Next, Cytoscape software (version 3.7.1) (http://www.cytoscape.org) was used to visualize the interaction relationships between genes in these hub nodes. Finally, a subnetwork was extracted from the PPI network using the Molecular Complex Detection (MCODE) application,[ a Cytoscape plugin that detected densely connected regions in networks that may represent molecular complexes. The parameters were set as follows. Degree cutoff: 2, node score cutoff: 0.2, cut style: haircut, k-core: 2, and max. Depth: 100. Genes in the PPI subnetwork were considered as hub genes for further analysis.

Hub gene identification and validation using TCGA data

The RNA-seq data and clinical information of PCa were downloaded from The Cancer Genome Atlas Project database (TCGA, https://cancergenome.nih.gov/). A total of 551 samples were included. The samples in the data were divided into different groups according to the clinical features of hub modules which we were interested in. Spearman correlation analysis was used to estimate the filtered gene expression in differential clinical traits. Then, event-free survival (EFS) estimate and statistical significance were performed using the Kaplan–Meier method, two-sided log-rank test, and univariate Cox regression analysis. A hazard ratio (HR) P value of <.05 was considered significant statistically. The Human Protein Atlas (http://www.proteinatlas.org) was applied for validating the immunohistochemistry (IHC) of latent hub genes.

Results

DEGs screening and WGCNA analysis

Totally 2688 DEGs (1539 up-regulated and 1149 down-regulated) were identified between PCa and matched normal adjacent tissue (Fig. 1A). The WGCNA analysis was performed on these DEGs. The tumor samples of GSE69223 were clustered using Pearson correlation method and the average linkage method (Fig. 1B). Meanwhile, a scale-free network was ensured with the power of β = 10. Totally 6 modules were identified via the average linkage hierarchical clustering (Fig. 1C). Both brown module (r = 0.58, P value = .02) and turquoise module (r = −0.61, P value = .02) were found to have the most conspicuous association with tumor grade (Fig. 1D). The expression levels of genes in the brown module were positively associated with tumor grade of PCa, while those of the turquoise module had the contrary relationship. These 2 modules were selected as hub modules for subsequent analysis.

Figure 1

DEGs screening and WGCNA analysis. (A) The volcano picture for DEGs. Green dots represented the downregulated genes. Black dots represented genes that were not differentially expressed, and the red dots represented the upregulated genes. (B) Clustering dendrogram of 15 tumor samples as well as the clinical features. The color intensity was proportional to higher tumor grade, higher tumor stage, and older age. (C) Dendrogram of all DEGs clustered based on a dissimilarity measure (1-TOM). (D) Heatmap of the correlation values between MEs and different clinical features of PCa (tumor stage, tumor grade, and age). Red for positive correlation and Blue for the negative correlation. P values were printed below the correlations. DEGs, differentially expressed genes. MEs = module eigengenes, PCa = prostate cancer, WGCNA = weighted gene co-expression network analysis.

Gene ontology and pathway enrichment analysis

With the similarly prominent connection of tumor grade, we combined the genes in these 2 modules for enrichment analysis. The genes were categorized into 3 functional groups (BP, CC, and MF). The results showed that DEGs in the BP group was mainly involved in oxidation–reduction process, cell adhesion, extracellular matrix organization, Wnt signaling pathway, and cytoskeleton organization; DEGs in the CC group were enriched primarily on cytoplasm, extracellular space, extracellular exosome, membrane especially plasma membrane; DEGs in the MF group were significantly improved in protein homodimerization activity, GTP binding, structural molecule activity, GTPase activity, and transferase activity. According to KEGG enrichment analysis, our findings indicated that DEGs mainly took part in focal adhesion, proteoglycans in cancer, vascular smooth muscle contraction, cGMP−PKG signaling pathway and ECM−receptor interaction (Fig. S3).

Based the cut-off criteria (showed in Fig. S4. Brown module:|MM| > 0.6 and |GS| > 0.2;Turquoise module: |MM| > 0.8 and |GS| > 0.2), totally 570 genes (286 genes in the brown module and 284 genes in the turquoise module) with high connectivity were identified as candidate hub genes. The number of highly correlated genes between the 2 modules would be so far apart that it might affect the eventual screening of hub genes if we used |MM| > 0.8 as cut-off criteria for brown module. Therefore, we used a lower cutoff for it. Figure 2 showed the 2 PPI networks created by importing these candidate hub genes of 2 hub modules into the STRING database, respectively. The 13 highest-scoring nodes were screened out by MCODE application: CCNB1, TTK, DBF4, WEE1, MCM3, MCM6, ACTG2, CNN1, TAGLN, MYL9, ACTA2, LMOD1, and FLNA (6 genes in the brown module and 7 genes in the turquoise module). These 13 genes were considered as hub genes.

Figure 2

Visualization of PPI networks with candidate hub genes in the brown and turquoise modules, respectively (286 genes in the brown module and 284 genes in the turquoise module). PPI = protein-protein interaction.

Based on the samples in TCGA data which were divided into 5 groups according to the tumor grades, Spearman correlation analysis showed that 2 up-regulated genes (CCNB1, TTK) and 2 down-regulated genes (ACTG2, CNN1) in PCa had the most significant correlation coefficient (|cor|>0.25), which meant there was a highly consistent variation tendency between the expression levels of these genes and the grades of PCa (Fig. 3). Meanwhile, Kaplan–Meier curves and the two-sided log-rank test revealed that lower expression of CCNB1, TTK, and higher expression of ACTG2, CNN1 were associated with better EFS (Fig. 4). Besides, the univariate Cox proportional hazards regression analysis revealed that these 4 genes were also correlative with EFS of PCa (Table 1). Moreover, based on IHC data, the protein levels of CCNB1 and TTK were slightly higher in tumor tissues than those in normal tissues, while that of ACTG2 in tumor tissues was significantly lower than those in the control group (Fig. S5).

Figure 3

The expression levels of (A) CCNB1, (B) TTK, (C) CNN1, and (D) ACTG2 in different tumor grades of PCa. PCa, prostate cancer.

Figure 4

Event-free survival curve of the four hub genes in PCa based on Kaplan–Meier analysis and log-rank test. Patients were divided into the high expression level group and the low expression level group based on quartile cutoff method (cutoff-high:75% and cutoff-low:25%). (A) CCNB1. (B) TTK. (C) CNN1. (D) ACTG2. PCa, prostate cancer.

Table 1

Univariate Cox analysis of the 4 hub genes using TCGA dataset.

Discussion

PCa is one of the most common malignant tumors in men, which seriously harms men's health. Current diagnostic methods such as PSA testing, digital rectal examination, and prostate biopsy are defective and might be rising in a number of non-malignant clinical situations. Thus, it is vital to search for novel non-invasive tumor markers with high sensitivity and high specificity. Many studies have shown that the development and application of multi-omics-based approaches including genomics, proteomics, and metabolomics in malignant tumors provide the possibility for this idea, which is conducive to improving the accuracy of early diagnosis and prognosis evaluation for PCa and reducing the number of unnecessary invasive examinations.[ For instance, the PCA3 (also called DD3) gene was reported remarkably increased for PCa tissue compared with the normal control, and PCA3 score of 20 could be regarded as the cutoff of avoiding unnecessary biopsies.[ Besides, 1 research figured out that the levels of Periostin, Cathepsin D, Lysosome-associated membrane glycoprotein 2 and Adipocyte plasma membrane-associated protein were detected highly specific in PCa serum with respect to the normal group, which were proved by a fast protocol based on highly sensitive mass spectrometric analysis by Selected Reaction Monitoring and TiO2 enrichment.[ What is more, another research summarized that 2 major metabolomics characteristics of PCa were an increased expression of genes in the tricarboxylic acid cycle and the accumulation of metabolic intermediates, which might guide us in exploring novel therapeutic targets for PCa.[

It is worth noting that genomics, as a relatively mature omics technology, has shown unique advantages in the exploration of biomarkers for malignant tumors. In the study, we performed an integrated bioinformatics analysis based on WGCNA towards PCa patients gene expression profiling, finally highlighting some key hub genes which could act as candidate biomarkers associated with the tumor progression and prognosis of PCa. We also made use of the RNA-seq data and clinical traits of PCa from the TCGA database, as well as IHC data from the Human Protein Atlas database for validation.

Both GO and KEGG pathway analysis in the combination of identified modules highly enriched in cell adhesion. Interestingly, numerous studies figured out that changes of adhesion-related proteins were related to high Gleason score and metastasis, such as E-cadherin[ and CD44 standard (CD44s).[

Additionally, another research revealed that glucose-regulated protein 78 was up-regulated in PCa and served as a modulator of cell adhesion markers, which contributes to metastatic spread.[ GO terms are also enriched in the Wnt signaling pathway that is of great interest to the PCa community. Wnt family contains several secreted glycoproteins, which could combine with some Frizzled transmembrane receptors to activating one of 2 principal Wnt pathways (i.e., “non-canonical” or “canonical”).[ Several studies demonstrated that WNT5a might serve as a tumor-inhibiting factor because its expression could induce apoptosis in PCa cell lines[ and predict for increased biochemical recurrence-free survival.[ Similarly, the expression level of WNT11 seems to be upregulated by androgen deprivation treatment,[ which contribute to neuroendocrine differentiation in PCa.[

Furthermore, we screened out 4 hub genes (CCNB1, TTK, CNN1, and ACTG2) with high connectivity from the identified modules through PPI networks. The expression patterns of these genes were correlated with the tumor grade (WGCNA), and the gene expression status was associated with prognosis. Our findings suggested that the mentioned genes might play a crucial role in the occurrence and development of PCa.

CCNB1, also named cyclin B1, regulates the G2 to M phase transition in mitosis. Its expression and clinical value have been examined in a variety of cancers. For example, overexpression of CCNB1 played an essential role in cell proliferation and tumor growth of colorectal cancer, which suggested that CCNB1 may be an effective biomarker in colorectal cancer.[ What is more, another bioinformatics analysis revealed that CCNB1 could serve as a potential therapeutic target in hepatocellular carcinoma.[ As for PCa, according to one previous research, androgen receptor could increase the expression level of CCNB1 via Akt phosphorylation, which contributes to the progression of PCa regulated by Jagged1.[ Besides, the high ratio of cyclin A and B to the proliferation marker Ki67 was proved as the most potent predictor of time to recurrence of PCa.[ Higher levels of CCNB1 sensitized PCa cells to apoptosis initiated by chemotherapy, supporting that overexpression of CCNB1 in PCa cells may be an effective prognostic biomarker for chemotherapy.[ In the present study, our results suggested that a high CCNB1 expression level could act as a good predictor of worse EFS.

TTK, also called Monopolar spindle one kinase (MPS1), is another checkpoint in mitosis and correlative with cell proliferation and cell invasion. Currently, overexpression of TTK was confirmed to occur in various cancers that harbor copy number gains (CNG) of chromosome 9p involving PD-L1.[ Furthermore, TTK was overexpressed within the incredible lions share of esophageal squamous cell carcinoma tissues but barely expressed in ordinary tissues,[ while another study demonstrated that high expression status of TTK corresponded positively with tumor grade and negatively with prognosis in glioblastoma.[ Nowadays, there are few types of researches regarding the role of TTK in PCa. Our study discovered that overexpression of TTK was associated with higher Gleason grade and worse prognosis of PCa, which was the first study revealing the relationship between TTK expression level and clinical outcome of PCa, as far as we know.

Functionally, the binding of CNN1 contributed to the stabilization of the 3D structure of actin filaments, which effectively resulted in the decrease of filament turnover and suppression of cytoskeleton remodeling.[ Unlike CCNB1 and TTK, CNN1 may play a tumor-suppressive role in many kinds of malignancy tumorsr.[ Loss of CNN1 also led to tumor progression and anoikis resistance in high-grade serous carcinoma, which originated from the fallopian tube epithelium.[ Although the functional importance of CNN1 on many other cancers has been proved, few studies have performed the specific analysis of CNN1 on its tumor-suppressive function in PCa. Based on the present study, we detected that the expression of CNN1 in PCa was down-regulated and related with high tumor grade and terrible clinical outcome, supporting CNN1 may be a candidate biomarker for progression and prognosis of PCa.

Recently, aberrant expression of ACTG2 (gamma smooth muscle 2) has been described in several different types of malignancy tumors. According to 1 previous study, a lower expression level of ACTG2 was correlative with worse disease-specific survival in leiomyosarcoma.[ On the contrary, it was also reported that ACTG2 was up-regulated in hepatocellular carcinoma, involved as a promoter gene in cell migration and distant metastasis.[ Higher expression status of ACTG2 was also found in colon tissue compared with colon carcinoma.[ Its expressions and effects seem to vary in different kinds of cancers. Although the molecular mechanism by which ACTG2 functions in tumors is not very clear and needs further research, it is reasonable to believe that ACTG2 could be a novel prognostic biomarker and candidate therapeutic target for PCa due to the results of our study.

Nevertheless, the limitations should be acknowledged for our research. Firstly, it is required further studies in vivo and in vitro for illuminating the accurate molecular mechanisms that affected the tumorigenesis and the progression of PCa. Secondly, to making our findings more convincing, a larger number of clinical samples were needed, which we will concentrate on in our future research.

Taken together, we identified 4 potential molecular biomarkers (CCNB1, TTK, CNN1, and ACTG2) for prognosing outcome and target therapy of PCa by applying a comprehensive bioinformatics analysis based on WGCNA, but more detailed functional exploration and validation of the molecular mechanisms associated with these genes are indispensable. Our findings might furnish some meaningful insights into individual therapy of PCa patients.

Author contributions

Conceptualization: Xuan Chen, Jingyao Wang, Yongqing Lai.

Data curation: Xiqi Peng.

Formal analysis: Xuan Chen, Jingyao Wang, Chunduo Zhang.

Funding acquisition: Yongqing Lai.

Investigation: Xuan Chen, Jingyao Wang, Xiqi Peng, Chunduo Zhang.

Methodology: Xuan Chen, Jingyao Wang, Xiqi Peng, Kaihao Liu, Xingzhen Zeng.

Project administration: Xingzhen Zeng, Yongqing Lai.

Resources: Xingzhen Zeng, Yongqing Lai.

Software: Xuan Chen, Jingyao Wang, Chunduo Zhang.

Supervision: Xingzhen Zeng, Yongqing Lai.

Validation: Xuan Chen, Xiqi Peng, Kaihao Liu.

Visualization: Xuan Chen, Kaihao Liu.

Writing – original draft: Xuan Chen.

Writing – review & editing: Jingyao Wang, Yongqing Lai.