Impact Factor 6.244 | CiteScore 3.9
More on impact ›


Front. Oncol., 02 February 2021 |

Identification of Key Gene Signatures Associated With Bone Metastasis in Castration-Resistant Prostate Cancer Using Co-Expression Analysis

Zhongxiang Yu1†, Hanlin Zou2†, Huihao Wang3, Qi Li4* and Dong Yu5*
  • 1Department of Orthopaedics, Shuguang Hospital Affiliated to Shanghai Traditional Chinese Medical University, Shanghai, China
  • 2Department of Orthopedics, Putuo Hospital Affiliated to Shanghai Traditional Chinese Medical University, Shanghai, China
  • 3Shi’s Center of Orthopedics and Traumatology, Shuguang Hospital Affiliated to Shanghai Traditional Chinese Medical University, Shanghai, China
  • 4Department of Oncology, Shuguang Hospital Affiliated to Shanghai Traditional Chinese Medical University, Shanghai, China
  • 5Center for Translational Medicine, Second Military Medical University, Shanghai, China

About 80–90% of castration-resistant prostate cancer (CRPC) patients would develop bone metastasis. However, the molecular mechanisms of bone metastasis are still not clear. This study aimed to detect the differences between the tumor and normal samples in bone after metastatic colonization. Four transcriptional datasets (GSE32269, GSE101607, GSE29650, and GSE74685) were obtained from the GEO database. 1983 differentially expressed genes (DEGs) were first identified between tumor and normal marrow samples in GSE32269. Most of the top 10 up-regulated DEGs are related with prostate cancer, and the top 10 down-regulated DEGs are mainly related with bone development. Seven co-expression modules were then detected based on the 1469 DEGs shared by the four datasets. Three of them were found highly preserved among the four datasets. Enrichment analysis showed that the three modules were respectively enriched in Cell adhesion molecules (CAMs), Leukocyte transendothelial migration and cell cycle, which might play significantly important roles in the tumor development in bone marrow. Ten, 17, and 99 hub genes for each module were then identified. And four genes (C3AR1, IL10RA, LY86, and MS4A6A) were detect to be tightly related to progression of bone metastatic CRPC. ROC curve was plotted and AUC was calculated to distinguish tumor and normal bone marrow samples as well as bone and non-bone metastatic CRPCs. The present study identified key genes and modules involved in bone metastatic CRPCs, which may provide new insights and biomarkers for understanding of the molecular mechanisms of bone metastatic CRPC.


Prostate cancer (PCa) is one of the most common cancers and the tenth most common cause of cancers related mortality in men in China (1). The rankings rise first in men in the developed countries (2). Castration-resistant prostate cancer (CRPC) is an advanced form of prostate cancer by disease progression following surgical or pharmaceutical castration. This process is not inevitable, which is usually companied by poor prognosis and reduced survival time. To be known, CRPC patients are also at high risk of developing metastases. The common sites are bone, lymph nodes, liver, lungs and brain. Bone is the most prominent site for metastases. About 80–90% of CRPC patients develop bone metastases (3). Bone metastases could lead to the disorder of bone metabolism and induce skeletal related events (SREs), such as pathological fracture, spinal cord compression and hypercalcemia, which not only reduce survival time and life quality, but also increase burden of treatment (4).

However, the molecular mechanisms of bone metastases are still not clear. A widely accepted mechanism is the ‘seed and soil hypothesis’, which describes an interaction between circulating tumor cell and microenvironment of bone tissue (5). The communication between bone and cancer cells is believed to be critical for the development and progression of bone metastases (6, 7). Most of researches focus on dissecting the process of initiation to development of distant metastasis, such as cancer cells migrating through the endothelial cells to gain access to systemic circulation via the tortuous and leaky tumor vasculature and cell signaling aberrations (8, 9). A set of marked differences were identified between metastases and primary tumors and the subgroups of bone metastasis were also detected by transcriptome or proteome analysis (1012). In addition, David A. Quigley et al. explore the genomic hallmarks and structural variation in metastatic PC, including bone metastatic CRPCs (13). However, these researches do not pay more attention on the state of tumor cells after metastatic colonization and also do not explore the differences between the tumor cells and normal cells in bone. This study aimed to identify the differences between tumor and normal bone marrow samples through differential expression analysis and weighted gene co-expression analysis. The identified key genes and modules will provide new insights for understanding of the molecular mechanisms and clinical treatment for bone metastatic CRPC.


Data Collection and Preprocessing

Four expression profile datasets containing CRPC bone metastasis were downloaded from the GEO database ( Dataset GSE32269 was chosen for further analysis with 29 CRPC bone metastatic marrow samples and four normal bone marrow samples, which was used for bone cancer significantly expressed genes selection and correlated modules detection. The other three datasets GSE101607, GSE29650, and GSE74685 were kept with only CRPC bone metastatic samples, which was used to validate and screen the truly significant and preserved bone cancer related modules. Detailed information of datasets was shown in Table 1.


Table 1 Datasets of gene expression profiles.

Before the analysis, all the raw data were reprocessed. Probes were mapped to the gene symbols. Empty probes and probes mapping to multiple genes were both discarded according to each annotation platform. If there were multiple probes that mapped to the same gene symbol, their mean values were considered as the gene expression value. The reprocessed data was normalized by the limma (linear models for microarray data) package in R (14).

Identification of Differentially Expressed Genes

The eBayes analysis was used to detect the differentially expressed genes (DEGs) between metastatic bone marrow samples and normal marrow samples in GSE32269 using limma package (14). The adjusted P-value <0.05 and |log-fold change|>1 were set as the threshold for DEGs screening.

Enrichment Analysis

R package clusterProfiler (15) was used for the Enrichment analysis. False discovery rate (FDR) < 0.05 was set as the threshold for the identification of significant GO-Enrichment terms and Pathway-Enrichment terms.

WGCNA Analysis

The co-expression network analysis was performed using weighted gene co-expression network analysis (WGCNA) (16). First, the soft threshold for network construction was selected, which is the lowest power for which the scale-free topology fit index curve flattens out upon reaching a high value. Second, the function blockwiseModules was used for one-step network construction and module detection. The module eigengene (ME) of each module and the correlation between MEs was then calculated. Thirdly, module preservation was calculated between GSE32269 and the other three datasets using the function modulePreservation (17). The comparability of two datasets is assessed by correlating measures of average gene expression and overall connectivity of two datasets. The higher the correlations of these properties, the better chance you will have of finding similarities between the two datasets at subsequent stages of analysis. Fourthly, the key node (hub gene) was determined by high intramodule connectivity of genes. The cut-off criteria was set |cor.geneModuleMembership| > 0.8. According to the intramodule connectivity, the detected hub genes were visualized using VisANT software (18). Finally, the study (19) containing mRNA and clinical data of 444 metastatic CRPC samples was used to validate the hub genes and subjected to survival analysis. The database GEPIA2 containing TCGA datasets (20) and the database Oncomine containing cancer microarray datasets (21) were used to validate the expression levels of hub genes.


DEG Identification for CRPC Bone Metastatic Patients

In order to detect the transcriptomic differences between CRPC bone metastatic marrow samples and normal marrow samples, the dataset GSE32269 with 29 CRPC bone metastatic marrow samples and four normal bone marrow samples was selected and downloaded from GEO databases. DEGs were identified using the limma package. 1983 DEGs were screened with the threshold of |logFC|>1 and p.adjust<0.05, as shown in >>Figure 1A, which contains 825 up-regulated genes and 1158 down-regulated genes for bone metastatic marrow samples (see Supplementary Table 1). The top 10 significantly expressed genes are KLK3, KRT18, EFNA1, SLC396A, NKX3-1, PGLYRP1, MGAM, RHD, GFI1, FAR2, of which the first five are up-regulated and the following five are down-regulated. The expression profiles of these DEGs were showed as a heatmap in Figure 1B. Enrichment analysis was further conducted. The result was shown in Figures 1C, D. The most enriched GO terms are neutrophil and leukocyte-associated terms. The top5 pathway terms are Malaria, Leukocyte transendothelial migration, B cell receptor signaling pathway, phagosome and chemokine signaling pathway.


Figure 1 The volcano, heatmap, GO and KEGG enrichment results of differentially expressed genes (DEGs) between tumor and normal cells in bone. (A) The volcano plot for DEGs. Grey dots represent genes which are not differentially expressed, red dots represent the upregulated genes, and the blue dots represent the downregulated genes. (B) The heatmap for DEGs. (C) The annotation of gene ontology function of DEGs using GO enrichment analysis. (D) The annotation of pathway function of DEGs using KEGG enrichment analysis.

WGCNA Analysis

Since the four datasets come from different platforms, we should ensure that the four datasets are comparable. First, we need to limit the analysis to genes that expressed among the datasets. The intersection was taken among the DEGs of GSE32269 and the genes of other three datasets. 1469 genes were selected, and the corresponding expression profiles of these genes in four datasets were then prepared. Second, the comparability of GSE32269 and other dataset was assessed by measuring the average gene expression and overall connectivity between two datasets (Figure 2). It’s clear to see that the correlations are positive and the p-value are significant in all cases, which suggests that the datasets are comparable.


Figure 2 The correlations of average gene expression (A) and overall connectivity (B) between GSE32269 and other three datasets (GSE101607, GSE74685, GSE29650).

Prior to gene co-expression network detection, the analysis of network topology for various soft-thresholding powers was performed to obtain relative balanced scale independence and mean connectivity. As shown in Figure 3A, power seven was the lowest power for which the scale-free topology fit index reaches 0.85. Based on this power, seven modules were generated as shown in Figure 3B. The largest module was the turquoise module, which contained 585 genes, the smallest module was the black module containing 49 genes. Averagely, each module contained 183 genes.


Figure 3 Identification of modules in the dataset GSE32269. (A) Network topology of different soft-thresholding powers. The left panel displays the influence of soft-thresholding power (x-axis) on scale-free fit index (y-axis). The right panel shows the influence of soft-thresholding power (x-axis) on the mean connectivity (degree, y-axis). (B) Clustering dendrogram showing eight modules that contain a group of highly connected genes. Each designated color represents a certain gene module.

Enrichment analysis was further performed to detect biological significance of each module as listed in Supplementary Table 2. In the top 5 terms of each module, Yellow, Turquoise and Brown module were mainly enriched in neutrophil-associated GO terms, which were all related with leukocyte mediated immunity. Red module had no significantly enriched pathways. Yellow and brown module shared an enriched pathway term, named Osteoclast differentiation, which is related with bone development. It’s worth noting that turquoise enriched pathways contain a set of signaling pathways, such as B cell receptor signaling pathway, chemokine signaling pathway, NF-kappa B signaling pathway, Fc epsilon RI signaling pathway and hematopoietic cell lineage. These are reported to be related with tumorigenesis. In the yellow module enriched pathways, Cell adhesion molecules are related with cancer invasion and metastasis. The green module is enriched with cell cycle-associated pathways.

Module Validation Among the Other Three Datasets

In order to detect whether these modules are preserved between the other three datasets, module preservation statistics were calculated using the function modulePreservation. The preservation Z-summaries was showed in Figure 4. We set the threshold Z>10 to screen the highly preserved modules. 3, 4, and 5 modules are separately found to be preserved in the dataset GSE29650, GSE74685 and GSE101607. And three modules (green, yellow, and turquoise) are shared and highly preserved in the three datasets, which were chosen for subsequent analysis.


Figure 4 The preservation Zsummary values of eight GSE32269 modules in other datasets (GSE29650, GSE74685 and GSE101607). The black horizontal line is the threshold to define the highly preserved modules among the four datasets.

Identification and Validation of Hub Genes

10, 17 and 99 hub genes were separately identified in the three preserved modules (Green, Yellow, Turquoise). The corresponding networks of hub genes were showed in Figure 5. The study containing mRNA data followed clinical information of 160 bone metastatic CRPC samples and 284 non-bone metastatic CRPC samples were subjected to survival analysis and regression analysis. Four hub genes (C3AR1, IL10RA, LY86, and MS4A6A) were identified to significantly associated with the overall survival (Figure 6). The patients with lower expression of the genes had a longer survival. However, the four genes have significantly higher expression level in bone compared to other non-bone metastatic tissues as showed in Figure 7. In CRPC patients with metastases, the bone metastases have the worst median progression to non-bone tissues metastases (22).


Figure 5 The visualization of hub genes in green module (A), yellow module (B), and turquoise module (C).


Figure 6 Survival analysis of hub gens with statistical significance (pvalue<0.05) in the dataset derived from Abida W’s study. Orange lines represent high expression of the hub genes and blue lines represent low expression.


Figure 7 Box-plot of expression values (FPKM) of the four hub genes between bone and non-bone tissues in the metastatic CRPC patients derived from Abida W’s study.

In addition, ROC curve analysis was implemented to evaluate the capacity of the hub genes to distinguish bone and non-bone metastatic tissues. AUC values for the four genes were greater than 0.6 (Figure 8).


Figure 8 ROC analysis of four hub genes in the dataset derived from Abida W’s study. Receiver operating characteristic (ROC) curves and area under the curve (AUC) statistics is to evaluate the capacity of distinguishing bone and non-bone metastatic tissues in the metastatic CRPC patients.


Development of bone metastases is a key and usual event in the progression of CRPC, which could lead to disorders of bone metabolism and skeletal related events. The median survival form men with bone metastases CRPC is approximately 1.5–2 years. The purpose of this study was to dissect the expression profile differences between the established metastatic tumor and normal bone marrow samples and then identified some key gene signatures and modules based on co-expression network analysis. These results will be helpful to deeply understand the molecular mechanisms of bone metastases and also provide candidate biomarkers for the prognosis prediction of bone metastatic CRPC patients.

The screened DEGs are found to be mainly related with prostate cancer and bone development. For example, among the top 10 up-regulated genes, KLK3 and KLK2, are highly enriched in prostate cancer, which are taken as effective biomarkers for diagnose and prognostic monitoring of prostate cancer (23). GOLM1 (24), FOLH1B (25), STEAP1 (26) and PLPP1 (27) are also identified as a candidate biomarker for prostate cancer. AGR2 expresses strongly in prostate tissue and show increased expression in prostate cancer (28). In a word, the up-regulated genes are mainly related with the tumorigenesis of prostate cancer. As for the top10 down-regulated genes, all of them are identified to be overexpressed in whole blood according to GTEx (29) and take part in embryonic development of blood and bone according to LifeMap Discovery (30). Therefore, the down-regulation of these genes would have effects on the function of bone or bone marrow, which might be genetic causes of SKE. These results indicated that the colonization in bone of metastatic CRPC cells not only keep the expression features of prostate cancer, but also induce new expression variations associated with bone. In another way, these results suggest the tissue specificity of DEGs, and the reliability of our results. Therefore, it is important to further dissect the expression differences between the established tumor and the normal bone marrow samples.

After a series of bioinformatic analysis, four hub genes identified from the three highly preserved co-expression modules among the four datasets were found to be tightly associated with overall survival in bone metastases CRPC patients. At present, there are no direct evidences to verify the functions of the four genes in prostate cancer or bone metastatic CRPC, but a set of researches showed that these genes were involved in the tumorigenesis and tumor proliferation in other cancers. It was shown that C3AR1 was significantly correlated with the overall survival In glioblastoma, which showed a longer survival time in the patients with lower expression of C3AR1 (31). In a recent study, over-expression of C3AR1 was proved to promotes HL-60 cell migration and invasion in vitro experiment (32). In other words, down-expression might decrease the migration and invasion capacities of tumor cells. Moreover, C3a, which binds to an orphan G protein-coupled receptor encoded by C3AR1, was reported as an immune regulator in the tumor microenvironment and act as insidious propagators of tumor growth and progression (33). In this respect, the down-regulation of C3AR1 might inhibit the process of tumor growth and progression. Therefore, these may be the reasons why the patients with lower expression of C3AR1 had good prognosis. IL10RA encodes a receptor for interleukin 10, which can inhibit the synthesis of proinflammatory cytokines. In colorectal cancer, the expression of IL10RA is found to be higher in healthy tissue than in the CRC tissue and showed association with the proliferation index, confirming the importance of IL10RA in the pathogenesis of CRC (34). However, increased level of IL10RA in the study population was not linked with overall survival time. In diffuse large B-cell lymphoma, IL10 receptor is highly expressed and predicts worse survival (35). Functional experiment showed that IL10 receptor plays an important role in IL10-JAK-STAT signaling pathway. Blocking IL10R would interrupt the IL10 autostimulatory loop and lead to cell death through cell cycle arrest and introduction of apoptosis. LY86 encodes the lymphocyte antigen 86, which may cooperate with CD180 and TLR4 to mediate the innate immune response to bacterial lipopolysaccharide (LPS) and cytokine production. LY86 was identified as a novel biomarker for the prediction of osteosarcoma prognosis and therapeutic targets (36). Moreover, healthy hematopoietic stem progenitor cells (HSPCs) can be transformed genetically by leukemia macrovesicles to over express LSC specific genes, which contains LY86, LRG1 and PDE9A and so on (37). These suggests that LY86 might play an important role in the transformation of localized normal bone marrow cells to cancer cells. MS4A6A encodes a member of the membrane-spanning 4A gene family, which display unique expression patterns among hematopoietic cells and nonlymphoid tissues. GWAS researches showed that MS4A6A is associated with heel bone mineral density and Alzheimer’s disease (38). MS4A6A was reported to be highly expressed in putative Tumor-associated macrophages (TAMs) populations. Previous reports suggest that TAMs may show an immunosuppressive M2 signature, which promotes tumorigenesis by suppressing immune surveillance and inducing angiogenesis, rather that the activating M1-type signature (39). In addition, a recent study found that high expression of MS4A6A was associated with poor progression-free survival of ovarian cancer (40), which is consistent with the result of this study. Therefore, this gene might take an important role in the colonization of metastatic cancer cells in bone marrow and tumorigenesis of localized bone marrow cells. In above-mentioned studies, high expression of the four genes are all significantly associated with poor prognosis, which is consistent with the performances. These will serve as important references to explore the molecular mechanisms of the genes on bone metastatic CRPC.

In our results, the four genes were all down-regulated in tumor bone marrow samples compared to the normal samples, which was different from the performances in other tumors described above, However, the four genes present consistency trends as this study in the lung squamous cell carcinoma according to TCGA datasets (20). Some of the four genes were also lowly expressed in ACC, COAD or DLBC (Figure 9). We also made a search of the four genes in Oncomine database (21) with parameters (Analysis Type: Differential analysis, Cancer vs. Normal analysis, Prostate cancer vs. Normal analysis; Data Type: mRNA). The results showed that these genes have no differences in expression between tumor and normal samples in most of prostate cancer datasets as listed in the Figure 10, which is consistent with the result in the TCGA prostate cancer dataset.


Figure 9 Box-plot of expression values (TMP) of the four hub genes between tumor and normal samples derived from the TCGA datasets. ACC, Adrenocortical carcinoma; COAD, Colon adenocarcinoma; DLBC, Lymphoid Neoplasm Diffuse Large B-cell Lymphoma; LUSC, Lung squamous cell carcinoma; PRAD, Prostate adenocarcinoma. Differential analysis between tumor and normal group was conducted using one-way ANOVA method. *pvalue < 0.05.


Figure 10 Comparison of the four hub genes across 21 analyses of prostate cancer datasets. The colored box denotes the gene’s percentile rank for that analyses. The more saturated the color, the higher the percentile rank. Red denotes over-expression; blue denotes under-expression. The datasets were listed on the right.

At present, a growing number of researches focus on the communication between tumor cells and bone stroma (41). Existing discoveries show that a vicious cycle of molecular crosstalk between tumor cells and the bone metastatic niche often take place in osteolytic bone metastasis (42). Targeting the bone metastatic niche is also evolving into a promising avenue for the prevention of bone metastatic relapse, therapeutic resistance, and other aspects of cancer progression (4345). Therefore, it is meaningful and important to dissect the differences between tumor cells and bone metastatic niche at different level, including transcriptome, which will be essential to explore the molecular mechanisms or interaction underlying the bone metastases and new clinical practice. Based on this consideration, this study has creatively used public data to dissect the expression differences between established tumor and normal bone marrow samples derived CRPCs. The first screened DEGs were involved in prostate cancer and bone development. And the followed illustrated four hub genes are not only associated with overall survival of bone metastatic CRPC samples, but also be capable of distinguishing bone metastases and non-bone metastases. These findings would greatly provide new insights and biomarkers for understanding of the molecular mechanisms and clinical treatment for bone metastatic CRPC.

Data Availability Statement

Publicly available datasets were analyzed in this study, these can be found in here: the NCBI Gene Expression Omnibus (GSE32269, GSE101607, GSE29650 and GSE74685).

Author Contributions

ZY and DY designed the study and drafted the manuscript. ZY, HZ, and HW performed all the data analysis. HW helped the preparation of figures and tables. QL and DY contributed to the writing of the manuscript. All authors contributed to the article and approved the submitted version.


This work was supported by Shanghai University of Traditional Chinese Medicine research grants (18LK038). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Conflict of Interest

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


This manuscript has been released as a pre-print at Research Square (46).

Supplementary Material

The Supplementary Material for this article can be found online at:

Supplementary Table 1 | All the differentially expressed genes between metastatic tumor and normal bone marrow samples in GSE32269.

Supplementary Table 2 | GO and KEGG enrichment analysis results for each module.


CRPC, Castration-Resistant Prostate Cancer; DEG, Differentially Expressed Gene, WGCNA, Weighted Gene Correlation Network Analysis; ME, module eigengene.


1. Chen W, Zheng R, Baade PD, Zhang S, Zeng H, Bray F, et al. Cancer statistics in China, 2015. CA Cancer J Clin (2016) 66:115–32. doi: 10.3322/caac.21338

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Siegel RL, Miller KD, Jemal A. Cancer statistics, 2016. CA Cancer J Clin (2016) 66:7–30. doi: 10.3322/caac.21332

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Wirth M, Tammela T, Cicalese V, Gomez Veiga F, Delaere K, Miller K, et al. Prevention of bone metastases in patients with high-risk nonmetastatic prostate cancer treated with zoledronic acid: efficacy and safety results of the Zometa European Study (ZEUS). Eur Urol (2015) 67:482–91. doi: 10.1016/j.eururo.2014.02.014

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Sathiakumar N, Delzell E, Morrisey MA, Falkson C, Yong M, Chia V, et al. Mortality following bone metastasis and skeletal-related events among men with prostate cancer: a population-based analysis of US Medicare beneficiaries, 1999-2006. Prostate Cancer Prostatic Dis (2011) 14:177–83. doi: 10.1038/pcan.2011.7

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Pedersen EA, Shiozawa Y, Pienta KJ, Taichman RS. The prostate cancer bone marrow niche: more than just ‘fertile soil’. Asian J Androl Asian J Androl (2012) 14:423–7. doi: 10.1038/aja.2011.164

CrossRef Full Text | Google Scholar

6. Hiraga T. Bone metastasis: Interaction between cancer cells and bone microenvironment. J Oral Biosci (2019) 61:95–8. doi: 10.1016/j.job.2019.02.002

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Zheng Y, Zhou H, Dunstan CR, Sutherland RL, Seibel MJ. The role of the bone microenvironment in skeletal metastasis. J Bone Oncol (2013) 2:47–57. doi: 10.1016/j.jbo.2012.11.002

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Wan LL, Pantel K, Kang YB. Tumor metastasis: moving new biological insights into the clinic. Nat Med Nat Med (2013) 19:1450–64. doi: 10.1038/nm.3391

CrossRef Full Text | Google Scholar

9. Weis SM, Cheresh DA. Tumor angiogenesis: molecular pathways and therapeutic targets. Nat Med Nat Med (2011) 17:1359–70. doi: 10.1038/nm.2537

CrossRef Full Text | Google Scholar

10. Ylitalo EB, Thysell E, Jernberg E, Lundholm M, Crnalic S, Egevad L, et al. Subgroups of Castration-resistant Prostate Cancer Bone Metastases Defined Through an Inverse Relationship Between Androgen Receptor Activity and Immune Response. Eur Urol (2017) 71:776–87. doi: 10.1016/j.eururo.2016.07.033

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Iglesias-Gato D, Thysell E, Tyanova S, Crnalic S, Santos A, Lima TS, et al. The Proteome of Prostate Cancer Bone Metastasis Reveals Heterogeneity with Prognostic Implications. Clin Cancer Res Clin Cancer Res (2018) 24:5433–44. doi: 10.1158/1078-0432.CCR-18-1229

CrossRef Full Text | Google Scholar

12. Djusberg E, Jernberg E, Thysell E, Golovleva I, Lundberg P, Crnalic S, et al. High levels of the AR-V7 Splice Variant and Co-Amplification of the Golgi Protein Coding YIPF6 in AR Amplified Prostate Cancer Bone Metastases. Prostate (2017) 77:625–38. doi: 10.1002/pros.23307

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Quigley DA, Dang HX, Zhao SG, Lloyd P, Aggarwal R, Alumkal JJ, et al. Genomic Hallmarks and Structural Variation in Metastatic Prostate Cancer. Cell (2018) 175:889. doi: 10.1016/j.cell.2018.10.019

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res (2015) 43:e47. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS (2012) 16:284–7. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinf (2008) 9:559. doi: 10.1186/1471-2105-9-559

CrossRef Full Text | Google Scholar

17. Langfelder P, Luo R, Oldham MC, Horvath S. Is my network module preserved and reproducible? PloS Comput Biol (2011) 7:e1001057. doi: 10.1371/journal.pcbi.1001057

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Hu ZJ, Chang YC, Wang Y, Huang CL, Liu Y, Tian F, et al. VisANT 4.0: Integrative network platform to connect genes, drugs, diseases and therapies. Nucleic Acids Res Nucleic Acids Res (2013) 41:W225–31. doi: 10.1093/nar/gkt401

CrossRef Full Text | Google Scholar

19. Abida W, Cyrta J, Heller G, Prandi D, Armenia J, Coleman I, et al. Genomic correlates of clinical outcome in advanced prostate cancer. Proc Natl Acad Sci U.S.A. (2019) 116:11428–36. doi: 10.1073/pnas.1902651116

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Tang Z, Kang B, Li C, Chen T, Zhang Z. GEPIA2: an enhanced web server for large-scale expression profiling and interactive analysis. Nucleic Acids Res (2019) 47:W556–60. doi: 10.1093/nar/gkz430

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Rhodes DR, Yu J, Shanker K, Deshpande N, Varambally R, Ghosh D, et al. ONCOMINE: a cancer microarray database and integrated data-mining platform. Neoplasia (2004) 6:1–6. doi: 10.1016/s1476-5586(04)80047-2

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Koo KC, Park SU, Kim KH, Rha KH, Hong SJ, Yang SC, et al. Prognostic Impacts of Metastatic Site and Pain on Progression to Castrate Resistance and Mortality in Patients with Metastatic Prostate Cancer. Yonsei Med J (2015) 56:1206–12. doi: 10.3349/ymj.2015.56.5.1206

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Fagerberg L, Hallström BM, Oksvold P, Kampf C, Djureinovic D, Odeberg J, et al. Analysis of the human tissue-specific expression by genome-wide integration of transcriptomics and antibody-based proteomics. Mol Cell Proteomics (2014) 13:397–406. doi: 10.1074/mcp.M113.035600

CrossRef Full Text | Google Scholar

24. Varambally S, Laxman B, Mehra R, Cao Q, Dhanasekaran SM, Tomlins SA, et al. Golgi protein GOLM1 is a tissue and urine biomarker of prostate cancer. Neoplasia (2008) 10:1285–94. doi: 10.1593/neo.08922

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Zhang Z, Wu H, Zhou H, Gu Y, Bai Y, Yu S, et al. Identification of potential key genes and high-frequency mutant genes in prostate cancer by using RNA-Seq data. Oncol Lett (2018) 15:4550–6. doi: 10.3892/ol.2018.7846

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Gomes IM, Arinto P, Lopes C, Santos CR, Maia CJ. STEAP1 is overexpressed in prostate cancer and prostatic intraepithelial neoplasia lesions, and it is positively associated with Gleason score. Urol Oncol (2014) 32:53 e23–9. doi: 10.1016/j.urolonc.2013.08.028

CrossRef Full Text | Google Scholar

27. Uhlen M, Zhang C, Lee S, Sjostedt E, Fagerberg L, Bidkhori G, et al. A pathology atlas of the human cancer transcriptome. Science (2017) 357(6352). doi: 10.1126/science.aan2507

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Maresh EL, Mah V, Alavi M, Horvath S, Bagryanova L, Liebeskind ES, et al. Differential expression of anterior gradient gene AGR2 in prostate cancer. BMC Cancer (2010) 10:680. doi: 10.1186/1471-2407-10-680

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Consortium GT. The Genotype-Tissue Expression (GTEx) project. Nat Genet (2013) 45:580–5. doi: 10.1038/ng.2653

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Edgar R, Mazor Y, Rinon A, Blumenthal J, Golan Y, Buzhor E, et al. LifeMap Discovery: the embryonic development, stem cells, and regenerative medicine research portal. PloS One (2013) 8:e66629. doi: 10.1371/journal.pone.0066629

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Di J, Li S, Li D, Xue H, Yang D, Liu Y. Mining TCGA database for genes of prognostic value in glioblastoma microenvironment. Aging (Albany NY) (2018) 10:592–605. doi: 10.18632/aging.101415

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Yan Q, Li Z, Cen J-N, Chen S-N, Pan J, Hu S-Y. Over-expression of C3AR1 Promotes HL-60 Cell Migration and Invasion. Zhongguo Shi Yan Xue Ye Xue Za Zhi (2017) 25:1–7. doi: 10.7534/j.issn.1009-2137.2017.01.001

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Sayegh ET, Bloch O, Parsa AT. Complement anaphylatoxins as immune regulators in cancer. Cancer Med (2014) 3:747–58. doi: 10.1002/cam4.241

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Zadka Ł, Kulus MJ, Kurnol K, Piotrowska A, Glatzel-Plucińska N, Jurek T, et al. The expression of IL10RA in colorectal cancer and its correlation with the proliferation index and the clinical stage of the disease. Cytokine (2018) 110:116–25. doi: 10.1016/j.cyto.2018.04.030

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Béguelin W, Sawh S, Chambwe N, Chan FC, Jiang Y, Choo J-W, et al. IL10 receptor is a novel therapeutic target in DLBCLs. Leukemia (2015) 29:1684–94. doi: 10.1038/leu.2015.57

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Shi Y, He R, Zhuang Z, Ren J, Wang Z, Liu Y, et al. A risk signature-based on metastasis-associated genes to predict survival of patients with osteosarcoma. J Cell Biochem (2020) 121:3479–90. doi: 10.1002/jcb.29622

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Razmkhah F, Ghasemi S, Soleimani M, Amini Kafi-Abad S. LY86, LRG1 and PDE9A genes overexpression in umbilical cord blood hematopoietic stem progenitor cells by acute myeloid leukemia (M3) microvesicles. Exp Hematol Oncol (2019) 8:23. doi: 10.1186/s40164-019-0147-8

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Ma J, Zhang W, Tan L, Wang H-F, Wan Y, Sun F-R, et al. MS4A6A genotypes are associated with the atrophy rates of Alzheimer’s disease related brain structures. Oncotarget (2016) 7:58779–88. doi: 10.18632/oncotarget.9563

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Martinez FO, Gordon S, Locati M, Mantovani A. Transcriptional profiling of the human monocyte-to-macrophage differentiation and polarization: new molecules and patterns of gene expression. J Immunol (2006) 177:7303–11. doi: 10.4049/jimmunol.177.10.7303

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Pan X, Chen Y, Gao S. Four genes relevant to pathological grade and prognosis in ovarian cancer. Cancer Biomark (2020) 29:169–78. doi: 10.3233/CBM-191162

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Ren G, Esposito M, Kang Y. Bone metastasis and the metastatic niche. J Mol Med (2015) 93:1203–12. doi: 10.1007/s00109-015-1329-4

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Ell B, Kang Y. SnapShot: Bone Metastasis. Cell (2012) 151:690–690.e1. doi: 10.1016/j.cell.2012.10.005

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Boyerinas B, Zafrir M, Yesilkanal AE, Price TT, Hyjek EM, Sipkins DA. Adhesion to osteopontin in the bone marrow niche regulates lymphoblastic leukemia cell dormancy. Blood Blood (2013) 121:4821–31. doi: 10.1182/blood-2012-12-475483

CrossRef Full Text | Google Scholar

44. Sison EA, McIntyre E, Magoon D, Brown P. Dynamic chemotherapy-induced upregulation of CXCR4 expression: a mechanism of therapeutic resistance in pediatric AML. Mol Cancer Res (2013) 11:1004–16. doi: 10.1158/1541-7786.MCR-13-0114

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Park SI, Liao J, Berry JE, Li X, Koh AJ, Michalski ME, et al. Cyclophosphamide creates a receptive microenvironment for prostate cancer skeletal metastasis. Cancer Res Cancer Res (2012) 72:2522–32. doi: 10.1158/0008-5472.CAN-11-2928

CrossRef Full Text | Google Scholar

46. Yu Z, Zou H, Wang H, Li Q, Yu D. Candidate biomarkers and gene modules investigation for bone tumor samples derived from castration-resistant prostate cancer bone metastasis patients using WGCNA. (2019). doi: 10.21203/rs.2.17291/v1

CrossRef Full Text | Google Scholar

Keywords: bone metastatic CRPC, differentially expressed genes, weighted gene co-expression network analysis, module, hub genes

Citation: Yu Z, Zou H, Wang H, Li Q and Yu D (2021) Identification of Key Gene Signatures Associated With Bone Metastasis in Castration-Resistant Prostate Cancer Using Co-Expression Analysis. Front. Oncol. 10:571524. doi: 10.3389/fonc.2020.571524

Received: 11 June 2020; Accepted: 14 December 2020;
Published: 02 February 2021.

Edited by:

Giuseppe Giaccone, Cornell University, United States

Reviewed by:

Paul B. Fisher, Virginia Commonwealth University, United States
Marco Rossi, University of Catanzaro, Italy

Copyright © 2021 Yu, Zou, Wang, Li and Yu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Dong Yu,; Qi Li,

These authors have contributed equally to this work