Comparative Analysis of Prostate Cancer Gene Regulatory Networks via Hub Type Variation.

BACKGROUND: Prostate cancer is one of the most widespread cancers in men and is fundamentally a genetic disease. Identifying regulators in cancer using novel systems biology approaches will potentially lead to new insight into this disease. It was sought to address this by inferring gene regulatory networks (GRNs). Moreover, dynamical analysis of GRNs can explain how regulators change among different conditions, such as cancer subtypes.
METHODS: In our approach, independent gene regulatory networks from each prostate state were reconstructed using one of the current state-of-art reverse engineering approaches. Next, crucial genes involved in this cancer were highlighted by analyzing each network individually and also in comparison with each other.
RESULTS: In this paper, a novel network-based approach was introduced to find critical transcription factors involved in prostate cancer. The results led to detection of 38 essential transcription factors based on hub type variation. Additionally, experimental evidence was found for 29 of them as well as 9 new transcription factors.
CONCLUSION: The results showed that dynamical analysis of biological networks may provide useful information to gain better understanding of the cell.

Introduction

The complexity and multigenic nature of cancer has necessitated various genomewide studies to achieve a systems-level understanding of the key genetic mediators involved in prostate cancer ( Most diseases are due to the collapse of cellular processes together with interaction networks (). Therefore, exploring the biological network for complex diseases provides an understanding of the functional alterations in chronic diseases (.

Network-based approaches contain many clinical applications to explore human diseases systematically. A better understanding of the effects of cellular networks on disease progression may lead to the identification of disease genes which, in turn, may offer better targets for drug development (. One focal point in cancer analysis is the reconstruction of Gene Regulatory Networks (GRN) (. However, cancer progression is a dynamic process with multiple stages; so, reconstruction of one static GRN may not be informative enough for the inference. Instead, reconstructing stage-specific GRNs during cancer progression and then comparing these GRNs would be beneficial to characterize the main genes and interactions involved in cancer progression.

The availability of genomewide gene expression data has helped develop various state-of-art GRN reconstruction methods (. These methods seek to identify putative gene regulatory interactions by assuming that alterations in the expression level of a regulator (such as a transcription factor) have a direct effect on the cognate regulated genes.

Empirical evidence of extensive GRN rewiring during cancer progression, along with the availability of GRN reverse engineering approaches, have inspired us to conduct a systematic investigation to characterize the topological changes that occur in a prostate cell's GRN during cancer progression. Therefore, in this study, an attempt was made to reveal candidate disease-associated genes and biomarkers for prostate cancer progression by integrative gene expression profiling and network analysis at a systematic level. In this way, four stage-specific GRNs were reconstructed based on a comprehensive prostate cancer gene expression dataset containing 171 different samples monitoring gene expression at different disease phases. Topological comparison of these four GRNs based on hub type variations recapitulates the previous findings about extensive GRN rewiring (. Through sub-network analysis, it is possible to identify significant genes which were supposed to change their hub type with highly relevant to specific phases of prostate cancer.

Enormous efforts have been made to identify biomarkers for various cancers by the analysis of different transcriptome data (. Moreover, there were similar studies for analysis of sub-networks or hub genes which had been helpful for the understanding of the metastasis of cancer at the molecular level (. Nonetheless, there are still few studies on identification of prostate cancer biomarkers for disease progression (. Therefore, in this study, a new integrative network-based approach was developed to detect party hubs and date hubs based on Degree and BN algorithms during cancer progression.

Our analysis led to identification of 38 important genes putatively involved in prostate cancer. Through extensive literature search, experimental evidences revealed the role of 76.3% of candidate genes in prostate cancer (Table 1). This level of experimental confirmation reflects the high accuracy of the proposed approach.

Table 1

The function of 29 critical transcription factors putatively involved in prostate cancer

Gene name	Function
AR	Androgen receptor gene transcripts are over-expressed in most metastatic prostate cancers (40).
ATF6	ATF6-mediated apoptosis is reported in many cancers such as prostate cancer (49, 50).
SMAD3	SMAD3 is an essential mediator of tumor suppression and apoptosis in prostate cancer (41).
GATA3	Prostatic GATA3 is involved in androgen regulation of the prostate-specific antigen gene (51).
HLF	HLF is involved in prostate cancer promotion and progression (52).
PBX1	The Pbx1-HoxC8 heterocomplex causes androgen-independent growth in prostate cancer (53).
FOXF1	FOXF1 has high expression in normal prostate and its expression decreases in prostate cancer (39).
TCF21	TCF21 methylation levels accurately discriminate bladder and prostate cancerous tissues from their normal counterparts (54).
STAT1	Progressive dysregulation of STAT1 in prostate cancer cells contributes to prostate tumor growth (55).
EGR3	EGR3 is significantly over-expressed in prostate cancer and is a candidate prognostic marker of poor outcome prostate cancer (56).
FOS	In human prostate cancer, up-regulation of FOS protein occurs in advanced diseases and correlates with MAPK pathway activation (57).
KLF6	KLF6 controls cell cycle progression and apoptosis and is usually inactivate in many cancers such as prostate, ovary and colon (37, 38).
FOXA1	High-level FOXA1 expression is associated with the development of metastatic prostate cancer and could be used to classify patients who are at higher risk for metastases (58).
ELK4	ELK4 plays important roles in cell growth regulation of prostate cancer cells. The level of the transcript correlates with the progression of the disease (59).
HOXC6	Down-regulation of HOXC6 due to decreased proliferation rates of cell line and the over-expression of it rescues the cells from apoptosis in prostate cancer (60).
VDR	Studies discriminated the impact of VDR ligands upon prostate cancer cell proliferation, differentiation, and apoptosis (42).
RARB	RARB gene methylation in prostate samples is associated with an increased risk of subsequent prostate cancer (61).
EZH2	Over-expression of EZH2 causes invasion and growth of prostate cells. It is also a good biomarker for detection of the problem at an advanced stage (62).
TFAP2A (AP-2)	Cytoplasmic expression of AP–2 is reduced in prostate cancer cells (63).
JUNB	JUNB has an important role in controlling prostate cancer and can be a target for cancer therapy (64).
SNAI2	Down-regulation of SNAI2 is associated with primary prostate cancers and is a negative regulator of proliferation in the cancer cells (65).
ZEB1	Cancerous phenotype in prostate cancer cells is associated with increased expression of ZEB1 (66).
MXD4 (MAD4)	MXD4 which is known to have antitumor properties is significantly up-regulated in treated PC (67).
MAZ	MAZ expression deregulation relates to progression of many cancer types and plays an important role in PCa pathogenesis (68).
HOXB13	Recurrent mutation in HOXB13 associates with an increased risk of hereditary prostate cancer (69).
SIM2	Studies suggested an involvement of SIM2 in prostate tumor cell and cancer progression (70).
INSM1	Investigation showed that INSM1 remarkably up-regulates at the advanced PC stages (71).
PLAGL1	PLAGL1 is a tumor suppressor gene that inhibits growth of tumor cells by controlling apoptosis and cell-cycle progression in prostate cancer (72).
FOXC1	It was indicated that FOXC1 links to androgen-associated growth status of prostate cancer (39).

Our study hereby demonstrates a useful approach for analysis of prostate cancer at the systematic level. For the genomewide investigations, this will be a fundamental attempt for future development of the translational medical informatics, which lead to better patient diagnostics with high-throughput data through systems biology (.

Materials and Methods

Network reverse engineering approaches

Reverse engineering of GRNs from whole genome data entails deciphering the underlying gene regulatory circuits by observing changes in gene expression profiles (. With advances in high-throughput technologies, several computational reverse engineering approaches using different statistical measures have been developed (, including information-theoretic network inference methods, which identify connections between genes by approximating the quantity of information common to any pair of genes. In the Dialogue on Reverse Engineering Assessment and Methods 5 (DREA-M5) challenge, the context likelihood relevance (CLR) algorithm by Faith et al ( had the best performance among information theory based approaches (

Briefly, CLR determines an interaction between two genes to be significant by estimating the significance of their Mutual Information (MI) value against a background distribution of the MI values of every other pair involving one of the two genes of interest. In this way, the significance level is dynamically determined for each interacting pair according to their expression profiles. Given a gene expression dataset and the significance scores calculated by CLR algorithm, the corresponding empirical False Discovery Rate (FDR) can be estimated by running the algorithm on randomly shuffled datasets. In this study, all GRNs were reconstructed using CLR with an FDR threshold of 0.05. It is important to note that CLR relies solely on the dependency between expression profiles to detect interactions. Consequently, the resulting network is a coexpression network; so, the GRN is extracted from this network by considering only interactions where at least one transcription factor is involved.

Prostate cancer microarray data

Prostate cancer microarray data were downloaded from the Gene Expression Omnibus (GEO) database, accession number GSE6919 (. This dataset contains 171 samples, including samples from normal prostate tissue free of any pathology (Normal with 18 samples), normal prostate tissue adjacent to tumors (Adjacent with 63 samples), primary prostate tumor tissue (Tumor with 65 samples), and metastatic prostate cancer (Metastasis with 25 samples). Microarray data were preprocessed and analyzed using the LIMMA package in R ( which was originally developed for differential expression analysis of microarray data. Quantile normalization and a moderated tstatistic were used to find differentially expressed genes. More detailed descriptions of the methods can be found in the original publications.

Network topological analysis

To predict the key regulators in the prostate cancer based on hub type variations, the stage-specific GRNs of prostate cancer were searched for transcription factors which either had a high number of connections or were bottleneck (. The bottleneck genes are important because if they are removed from a network, the network will be disrupted, as they are major intersections between clusters in the network (. To find such genes, all constructed GRNs were topologically analyzed following the same rules using cyto-Hubba package ( which is a plugin for one of the most useful structural analysis software, Cytoscape (. Cyto-Hubba is used to detect the critical nodes of biological networks with many topological algorithms such as Degree, Bottleneck (BN), Maximum Neighborhood Component (MNC), Density of Maximum Neighborhood Component (DMNC), and a Double Screening Scheme (DSS) (. As mentioned before, Degree and BN were used to find the top ranked genes in all GRNs based on hub type variation which was the ultimate goal in this research.

Results

Stage specific network reconstruction of prostate cancer

Using the CLR algorithm, four independent networks related to the four different cell stages were reconstructed (Normal, Adjacent, Tumor, and Metastasis). The metastasis GRN had the lowest number of interactions with 2505 interactions while the other three GRNs had around 3000 interactions each. Additionally, topological analysis of the GRNs revealed that all four networks exhibited the small-word property ( and scales-free architecture ( which are the well-known characteristics of most biological networks (Figure 1). All four reconstructed GRNs were mainly composed of the same set of genes; however, the conserved interactions among these four networks were very low and the metastasis network had the most unique interactions (Figure 2).

Figure 1

The architecture of gene regulatory networks. All four networks (normal, adjacent, tumor and metastasis networks) follow the well-known characteristics of most biological networks; A) scale-free architecture defined as few highly connected genes (hubs) that link the other less connected genes to the network; B) small-word property which means any two genes in the network can be connected by relatively short paths through all interactions.

Figure 2

Number of interactions. This figure shows the GRNs for 14 TFs (orange nodes) that change their interaction numbers dramatically during cancer progression; A) Normal stage; B) Adjacent stage; C) Tumor stage; D) Metastasis stage which reflects the high level of rewiring of gene regulatory interactions.

Detection of essential transcription factors involved in the prostate cancer

Considering the importance of hub and bottleneck proteins in the structure of GRNs, the 50 highest-ranked genes were identified for each stage-specific GRN based on their degree and bottleneck scores, separately. Top 50 genes were selected based on previous studies that showed the highest percentage of critical proteins found in top 50 ones based on Degree and BN algorithm (. Although there were four conditions that resulted in detection of 200 genes, 144 of 200 genes had overlap during various conditions. Consequently, 56 unique candidate genes were selected for further analysis.

In each GRN, these 56 genes were categorized based on their degree and bottleneck scores in four groups: 1) Hub-NonBottleneck: genes with high degrees and low bottleneck scores are putative party hubs (; 2) Hub-Bottleneck: genes with high degrees and high bottleneck scores are putative date hubs (; 3) NonHub-Bottleneck: genes with low degrees and high bottleneck scores; 4) NonHub-NonBottleneck: genes with low degrees and low bottleneck scores.

The results showed hub type variation for 38 genes across different stages, whereas 18 other genes were functionally conserved as date hubs under all conditions (Table 2).

Table 2

56 transcription factors showed different topological characteristics in different stages

Gene name	N	A	T	M	Gene Name	N	A	T	M
HOXB13	PH	DH	PH	DH	ZNF516	DH	DH	DH	DH
ZEB1	PH	DH	DH	DH	GATA3	DH	DH	DH	NB
SIM2	PH	DH	DH	NB	ZNF423	DH	DH	DH	DH
TFAP2A	PH	DH	DH	DH	KLF6	DH	NN	NN	DH
ZNF205	PH	DH	DH	DH	ZNF91	DH	DH	DH	DH
EZH2	PH	NN	DH	PH	ETS2	DH	DH	DH	DH
MXD4	PH	DH	DH	DH	NFYA	DH	DH	DH	DH
ZNF146	PH	DH	DH	DH	MEIS1	DH	DH	DH	DH
CAMTA1	PH	NB	NN	DH	SMAD3	DH	NB	DH	DH
MAZ	PH	DH	DH	DH	MEIS2	DH	DH	DH	DH
PLAGL1	PH	DH	DH	PH	STAT5A	DH	DH	DH	DH
HOXC6	DH	PH	DH	PH	PRRX1	DH	DH	DH	DH
FOS	DH	PH	DH	DH	NHLH2	DH	NN	DH	DH
ELK4	DH	PH	NN	DH	FOXN3	DH	DH	DH	DH
NKX2-2	DH	DH	PH	DH	EGR3	DH	DH	DH	NN
MNX1	NB	NN	PH	DH	ID1	DH	DH	DH	DH
AR	DH	DH	DH	PH	NR4A1	DH	DH	DH	DH
PBX1	DH	DH	DH	PH	NR1H2	DH	DH	DH	DH
FOXA1	DH	DH	DH	PH	FOXF1	DH	DH	DH	NN
INSM1	NB	DH	NN	PH	ATF6	DH	DH	NB	DH
MEF2C	DH	DH	DH	DH	STAT6	DH	DH	DH	DH
STAT1	DH	DH	DH	NB	VDR	DH	DH	NN	DH
NR3C1	DH	DH	DH	DH	FOXC1	NB	DH	DH	DH
HLF	DH	DH	DH	NB	STAT2	NB	DH	NB	DH
TP63	DH	DH	DH	DH	JUNB	NN	DH	DH	DH
TCF21	DH	DH	DH	NB	RARB	NB	DH	DH	NB
ISL1	DH	NB	DH	NB	SNAI2	NB	NN	DH	NB
EGR1	DH	DH	DH	DH	ZNF529	NN	NN	NB	DH

N) Normal; A) Adjacent; T) Tumor; M) Metastasis; DH) Date Hub; PH) Party Hub; NB) Nonhub-Bottleneck; NN) Nonhub-Nonbottleneck.

Sub-networks consisting of the first neighborhoods of the 38 critical bottleneck transcription factors were extracted and compared, revealing changes in the number of interactions and gene targets between the stages. For some transcription factors such as STAT1, AR, HLF, ZEB1, TCF21, ISL1, KLF6, HOXB13, SIM2 and FOXA1, the number of interactions in the metastasis stage decreased dramatically (Table 3); the interaction numbers of other transcription factors, such as ZNF-529, FOXC1, MNX1, and JUNB, increased considerably in the metastasis stage, indicating network rewiring (Table 3). These 14 transcription factors (Figure 2) showed dramatic changes in their number of interactions (fold change ≥2) during the cancer progression (Table 3).

Table 3

The number of interactions for 14 out of 38 transcription factors changed dramatically from normal to metastasis stage (fold change >2)

Gene name	Normal	Adjacent	Tumor	Metastasis	Fold change
HOXB13	143	176	43	52	2.75
ZEB1	135	137	74	33	4.09
SIM2	120	132	115	25	4.80
TFAP2A	60	30	25	63	0.95
ZNF205	56	38	28	44	1.27
EZH2	46	24	81	81	1.76
MXD4	43	31	32	69	1.60
ZNF146	42	45	60	72	1.71
CAMTA1	41	22	21	37	1.11
MAZ	39	64	56	46	0.85
PLAGL1	26	36	39	30	0.87
HOXC6	34	53	78	47	1.38
FOS	63	49	69	65	0.97
ELK4	35	34	17	33	0.94
NKX2-2	28	33	43	30	0.93
MNX1	26	20	29	70	2.69
AR	103	61	35	48	2.15
PBX1	45	64	117	37	0.82
FOXA1	108	153	41	34	3.18
INSM1	21	26	23	32	0.66
STAT1	113	97	122	16	7.06
HLF	123	120	86	25	4.92
TCF21	59	102	58	25	2.36
ISL1	65	26	68	22	2.95
GATA3	42	34	68	22	1.91
KLF6	76	26	20	37	2.05
SMAD3	41	18	28	40	1.03
NHLH2	27	18	32	34	0.79
EGR3	32	35	40	26	1.23
FOXF1	32	151	129	26	0.81
ATF6	44	27	17	40	1.1
VDR	30	27	18	33	0.91
FOXC1	14	42	68	33	2.36
STAT2	23	32	19	34	0.68
JUNB	21	50	47	48	2.29
RARB	21	31	27	26	0.81
SNAI2	25	25	34	24	1.04
ZNF529	10	26	24	42	4.2

Discussion

To reconstruct cell stage specific GRNs, an attempt was made to focus on the available comprehensive transcriptome dataset, originally published (. This dataset was generated by sampling from four different types of prostate tissues including normal cells (Nor-mal), normal cells adjacent to cancer cells (Adjacent), primary tumor cells (Tumor), and metastatic cells (Metastasis). In our approach, genes were analyzed and prioritized based on the transcriptome data. Hence, it was possible to make reliable predictions only for genes with altered expression level across various conditions. To focus on these genes, only up/down-regulated genes were considered (fold change ≥2 and p < 0.05) in at least one cell stage (978 genes). Also, enrichment of known cancer genes was checked among this set by using a previously curated list of 555 high confidence cancer genes, originally published (. 100 genes were collected and identified as mediators in metastatic prostate cancer from (, and 276 genes were added and annotated as either a cancer pathway or prostate cancer gene in the KEGG database. It was found that the cancerrelated genes were about 1.55-fold (hypergeometric two-tailed test, p = 2.52E-6) and the prostate cancerrelated genes were about 2.22-fold enriched (hypergeometric two-tailed test, p = 6.07E-6) in our selected gene set which were fluctuated during prostate cancer.

To identify master regulators and their associated circuits governing cell-specific behavior in each state, the GRNs of prostate cells were compared in different stages with each other. Because the CLR algorithm merely relies on the similarity of expression patterns to infer interactions, constructed networks in this step contain both regulatory interactions (interactions between regulated genes and their putative regulators) as well as interactions between co-regulated genes (non-regulatory interactions). Hereafter, these networks are called coexpression networks. To extract gene regulatory interactions from these networks, only interactions involved at least one human transcription factor were considered and a list of them were extracted (. These networks are referred to as GRNs.

To predict the key genes in the prostate cancer, an attempt was made to find the stage-specific coexpression networks of prostate cancer for high connectivity (hub) or bottleneck genes (. Hub and bottleneck properties are considered important centrality indices because they are major intersections between clusters in the network and if they are removed from a network, the network will be disrupted (. Han et al suggested the existence of two types of protein hubs in the protein-protein interaction networks, namely party hubs and date hubs (. Although both interact with many proteins, the difference is that party hubs are proteins that interact with many other proteins simultaneously, whereas date hubs interact with their partners asynchronously (. By definition, the bottleneck proteins are responsible for the interconnection of clusters in the network, and thus bottlenecks with high degrees are most likely to be date hubs which contain groups of genes that assist in presenting common functions (. The obtained results recapitulate previous findings in which some active sub-networks contained regulatory interactions were supplanted by new interactions which changed their degrees during different conditions (.

The result also reflected the high level of rewiring of gene regulatory circuits during cancer progression, as suggested elsewhere (. As shown in table 3, for example, more than 2-fold decrease in the number of interactions for KLF6 was observed which controlled cell cycle progression and apoptosis. Indeed, experimental data suggest that KLF6 is inactivated in many cancers such as prostate, ovary and colon (. On the other hand, consistent with more than 2-fold increase in the number of interactions (from 14 in normal stage to 33 in the metastasis stage) for FOXC1 (Table 3), it was indicated that this gene is linked to androgen dependent growth of prostate cancer (.

Our result led to identification of 38 transcription factors which were bottleneck and changed their interaction during cancer progression. Although the functional role of some famous transcription factors such as AR, SMAD3 and VDR are well known as genes linked to prostate cancer (, the 9 transcription factors (CAMTA1, ISL1, MNX1, NHLH2, NKX2-2, STAT2, ZNF146, ZNF205, and ZNF529) are new candidates that may have critical roles in prostate cancer based on topological significance and regulatory changes during cancer progression. Among the remaining 9 transcription factors, 5 of them were associated with other cancer types. ZNF146, CAMTA1, NKX2-2, MNX1, and ISL1 are most prominent in colorectal cancer, neuroblastoma, Ewing's sarcoma, leukemia and breast cancer, and bladder cancers, respectively (. No evidence could be found to show the relationship between the 4 remaining transcription factors and any type of cancer.

Conclusion

In this paper, an accurate network-based framework for the analysis of transcriptome data was presented. The analysis of prostate state specific GRNs revealed 38 transcription factors which are critically important for prostate cancer progression. Also, 14 transcription factors were identified to be linked putatively to prostate cancer metastasis stage, so they would be used as key factors for future research in the field of cancer studies. Additionally, experimental evidences revealed the role of 29 of candidate transcription factors in prostate cancer.

The low number of predictions and high degree of overlap with previously known events in the prostate cancer demonstrate the high efficiency of our approach. In addition, the low number of predicted gene sets makes it easy to design follow up experiments to validate the results. In this study, it is believed that the results may provide critical information to gain better understanding of networks dynamics in the cell through complex diseases such as cancer.