Functional analysis of gene expression profiling-based prediction in bladder cancer.

The present study aimed to analyze the modification of gene expression in bladder cancer (BC) by identifying significant differentially expressed genes (DEGs) and functionally assess them using bioinformatics analysis. To achieve this, two microarray datasets, GSE24152 (which included 10 fresh tumor tissue samples from urothelial bladder carcinoma patients and 7 benign mucosa samples from the bladder), and GSE42089 (which included 10 tissues samples from urothelial cell carcinoma patients and 8 tissues samples from the normal bladder), were downloaded from the Gene Expression Omnibus database for further analysis. Differentially expressed genes (DEGs) were screened between benign the mucosa and control groups in GSE24152 and GSE42089 datasets. Gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) analysis were performed on overlapping DEGs identified in GSE24152 and GSE42089. Protein-protein interaction (PPI) networks and sub-networks were then constructed to identify key genes and main pathways. GO terms analysis was also performed for the selected clusters. In total, 1,325 DEGs in GSE24152 and 647 DEGs in GSE42089 were screened, in which 619 common DEGs were identified. The DEGs were mainly enriched in pathways and GO terms associated with mitotic and chromosome assembly, including nucleosome assembly, spindle checkpoint and DNA replication. In the interaction network, progesterone receptor (PGR), MAF bZIP transcription factor G (MAFG), cell division cycle 6 (CDC6) and members of the minichromosome maintenance family (MCMs) were identified as key genes. Histones were also considered to be significant factors in BC. Nucleosome assembly and sequence-specific DNA binding were the most significant clustered GO terms. In conclusion, the DEGs, including PGR, MAFG, CDC6 and MCMs, and those encoding the core histone family were closely associated with the development of BC via pathways associated with mitotic and chromosome assembly.

Introduction

Bladder cancer (BC) is a heterogeneous disease with a variable disease history and is one of the most common genitourinary malignancies worldwide (1). A total of 72,570 new cases of urinary BC were diagnosed in the United States and 15,210 patients succumbed to BC in 2013 (2). The most common pathological type of BC is urothelial cell carcinoma, and the 5-year survival rate of BC is ~77% in the United States (3). BC is a disease with multifactorial etiology; radical cystectomy with pelvic lymphadenectomy is considered to be standard therapy (4), although the use of radiotherapy is considered as an alternative, particularly in more frail patients (5). However, these therapies are limited in their effectiveness, as BC has high recurrence and mortality rates, meaning that greater understanding of the course of BC development is urgently required.

A previous study suggested that different mechanisms have evolved to respond to specific phenotypic alterations in genes and cellular pathways involved in BC, and certain genetic variations in major tumor-associated pathways were proven to induce BC (6). Shen et al (6) assessed the differentially expressed genes (DEGs) and interacting pathways in BC using bioinformatics analysis, given that the genes encoding activator protein 1 (AP-1) and nuclear factor of activated T cells were key in BC. Zhou et al (7) analyzed the gene expression in human BC samples using the GSE42089 microarray dataset, identifying a set of genes associated with mitotic spindle checkpoint dysfunction as being key in BC. However, the fact that only one microarray dataset was used by each of the above studies may prove to be a limitation to the analysis described.

In the present study, the gene expression profiles examined had to use the same platform as GSE42089 (Affymetrix GeneChip Human Genome U133 Plus 2.0 Array) and had to be samples composed of bladder cancer specimens and normal bladder tissue. The only other dataset that fit these criteria was GSE24152. Therefore, two microarray profiles, GSE24152 and GSE42089, were used for integrated analysis of gene expression modification in BC in the present study; the use of the two gene expression profiles enabled a more reliable conclusion to be drawn. DEGs were identified, and gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed. Protein-protein interaction (PPI) networks and sub-networks were also constructed to identify the key genes and main pathways involved in BC. Using the aforementioned bioinformatics methods, the modification of gene expression in BC was analyzed by identifying significant DEGs and pathways, which may provide novel insight for the etiology and treatment research of BC.

Materials and methods

Microarray data

Two gene expression profiles, GSE24152 (8) and GSE42089 (7), were downloaded from the Gene Expression Omnibus database in National Center for Biotechnology Information (http://www.ncbi.nlm.nih.gov/geo/), based on the GPL6791 and GPL9828 platforms in the Affymetrix GeneChip Human Genome U133 Plus 2.0 Array, respectively. The GSE24152 dataset included 17 samples, of which 10 were fresh tumor tissue samples collected from patients with urothelial carcinoma of the bladder and 7 were benign mucosa samples from the bladder. The microarray GSE42089 dataset consisted of 18 samples, of which 10 samples were tissues from urothelial cell carcinoma and 8 were normal bladder tissue.

Data preprocessing and DEGs analysis

The robust multiple average algorithm in the affy package (9) was used to normalize the microarray data and box plots were then generated. The microarray data from GSE24152 and GSE42089 were divided into two groups, a bladder carcinoma group and a normal group. Using the Limma package (10), the probe-level data of two sets were converted into expression measures and DEGs were identified in the bladder carcinoma group compared with the control group in GSE24152 and GSE42089. A Venn diagram was generated using VennDiagram package (11) to screen common DEGs in GSE24152 and GSE42089 for further analysis. A false discovery rate <0.01 and a |log2FC (fold-change) |>0.5 was used as the threshold. Heat maps were generated using the heatmap.2 function in ggplot 2 (12) to display the relative expression differences of DEGs. In addition, the cor.test function was used to evaluate the correlation coefficient of two groups of DEGs.

GO and pathway enrichment analysis

GO terms analysis is a widely used approach for studying large-scale genomic or transcriptomic data in function that consists of three terms: Biological process, cellular component and molecular function categories (13). KEGG is a widely used collection of online databases that deals with genomes, enzymatic pathways, and biological chemicals (14). In the present study, GO analysis and KEGG pathway enrichment were performed on the three categories of the screened DEGs using the DAVID (15) analysis tool, with a statistically significant cut-off of P<0.01.

PPI network and sub-network construction

Cytoscape (16) is an open-source bioinformatics software platform used for the visualization of molecular interaction networks and integrating with gene expression profiles and other state data. In the present study, Biosgenet in Cytoscape (17) was used to predict and visualize the interactions between selected DEGs and proteins using the BIND database (18) with a statistically significant cut-off of P<0.05. Significant sub-networks (clusters) were extracted from the PPI network using ClusterOne in Cytoscape (19) with a statistically significant cut-off of P<0.01. GO terms analysis was also performed for the selected clusters.

Results

DEGs selection

The variations in raw expression data were normalized following preprocessing (Fig. 1). A total of 1,325 DEGs were screened from the GSE24152 dataset and 637 DEGs were screened from the GSE42089 dataset (Fig. 2A and B depicts volcano plots of the two microarrays). A total of 619 common DEGs were identified by a Venn diagram, including 313 upregulated genes and 306 downregulated genes. From the heat maps of the 619 common DEGs generated (Fig. 2C and D), the DEGs could significantly distinguish the bladder carcinoma dataset from normal dataset.

Figure 1.

Box plots of data normalization. (A) Box plot of GSE24152 data normalization. (B) Box plot of GSE42089 data normalization. The x-axis represents samples; the y-axis represents gene expression values. The whiskers of the box plots represent the sample median.

Figure 2.

Volcano plots and heat maps of DEGs of two microarray profiles. Volcano plots of DEGs in (A) GSE24152 and (B) GSE42089. Heat map of DEGs in (C) GSE24152 and (D) GSE42089. DEG, differently expressed gene.

GO and KEGG enrichment analysis of DEGs

GO and KEGG enrichment analysis were performed with a cut-off of P<0.01, obtaining 74 GO terms and 4 KEGG pathways were obtained. The main GO terms and KEGG pathways of 619 common DEGs are listed in Table I. According to the results, DEGs were mainly enriched in chromosomal assembly and cell cycle pathways, including DNA replication, and the main GO terms were also associated with mitosis, including mitotic sister chromatid segregation and spindle checkpoint.

Table I.

Functional enrichment of differentially expressed genes.

Category	Term	P-value	Fold enrichment
BP	Mitotic sister chromatid segregation	1.90×10−4	9.8
BP	Negative regulation of cell cycle process	1.90×10−3	11.0
BP	Spindle checkpoint	9.70×10−3	16.3
CC	Condensed chromosome, centromeric region	1.70×10−9	9.8
CC	Condensed chromosome kinetochore	2.90×10−9	10.5
CC	Midbody	1.20×10−3	13.5
MF	Chromatin binding	4.50×10−3	3.8
MF	Microtubule motor activity	5.30×10−3	5.2
MF	Damaged DNA binding	1.70×10−3	7.3
KEGG	Progesterone-mediated oocyte maturation	1.80×10−4	5.5
KEGG	Cell cycle	1.10×10−8	6.1
KEGG	DNA replication	1.70×10−4	9.1

BP, biological process; CC, cellular component; MF, molecular function; KEGG, Kyoto Encyclopedia of Genes and Genomes.

A PPI network was constructed for all DEGs (Fig. 3) and clustering analysis was performed based on this PPI network. In total, 6 sub-networks (clusters) were obtained (Table II and Fig. 4). Progesterone receptor (PGR) was a key gene in cluster 1 and the proteins enriched in cluster 1 were all histone proteins. In cluster 2, MAF bZIP transcription factor G (MAFG) was identified as a key gene and it was also detected in cluster 5. Cell division cycle 6 (CDC6) was a key gene in clusters 3 and 4. Minichromosome maintenance complex component (MCM) family members, including MCM2, MCM4, MCM7 and MCM10, were also key genes in clusters 3 and 4. Nucleosome assembly and sequence-specific DNA binding were significant GO terms in cluster 1 and cluster 2, respectively. No GO terms were obtained for clusters 3–6.

Figure 3.

Protein-protein interaction network of DEGs. Pink boxes depict DEGs; blue boxes depict target proteins; pink arrows depict interactions between genes and target proteins; blue arrows depict interactions between proteins. DEG, differently expressed gene.

Table II.

Sub-networks obtained by clustering analysis.

Cluster	Nodes	Density	Quality	P-value	DEGs	Significant GO term
1	29	0.608	0.946	<0.001	PGR, CHAF1B, ASF1B,	Nucleosome assembly
2	16	0.653	0.842	9.18×10−5	HLF, MAFG	Sequence-specific DNA binding
3	15	0.514	0.635	4.41×10−5	CDC6, MCM2, MCM7, MCM10	None found
4	13	0.513	0.494	2.00×10−3	CDC6, MCM2, MCM10, MCM7, MCM4	None found
5	12	0.515	0.576	3.00×10−3	MAFG	None found
6	4	0.500	1.000	1.00×10−2	EFNA4	None found

DEGs, differentially expressed genes; GO, gene ontology.

Figure 4.

Sub-networks of 6 clusters. Pink boxes depict differently expressed genes; blue boxes depict target proteins; pink arrows depict interactions between genes and target proteins; blue arrows depict interactions between proteins.

Discussion

BC is a common malignancy that requires a high degree of surveillance, owing to the frequency of recurrence and the poor clinical outcome of invasive disease (20). Bioinformatic analysis of BC cells at the gene level can provide novel insights into this disease. In the present study, using the microarray datasets GSE24152 and GSE42089, significant DEGs in BC were identified. In the significant extracted sub-networks, PGR, MAFG, CDC6 and members of the MCM family were identified as key genes in BC; histones were also considered to have major functions in BC. The main GO terms and pathways were those associated with the cell cycle and chromosome assembly, including nucleosome assembly, spindle checkpoint and DNA replication.

PGR encodes a member of the steroid receptor superfamily, which mediates the physiological effects of progesterone (21). In the present study, PGR was an important gene in a cluster that was regulated by a set of transcription factors, revealing the potential significance of PGR in BC. Men are more frequently affected by BC than women, indicating that hormones and their receptors may function as regulatory factors (6). Miyamoto et al (22) clarified that the androgen receptor (AR) was involved in BC. As AR and PGR are determinant modulators of gonadal sex hormones, PGR may perform similar functions in BC with AR, via the progesterone-mediated oocyte maturation pathway, which was enriched in this study (23).

Histones are the main structural proteins associated with DNA in eukaryotic cells. Histones are divided into two groups, core histones and the nucleosomal histones (24). Core histones are some of the most highly conserved proteins in eukaryotes and have key roles in the organization of DNA folding (25). The altered patterns of histone modifications in various human cancer types have been studied extensively in recent years (26). Schneider et al (27) demonstrated that global histone modification levels were lower in BC than in normal urinary tissue. The conserved histone H2A has been reported to be overexpressed in BC cells and contributes to the activation of cancer-associated transcription pathways (28). In the present study, a set of core histones was clustered in cluster 1, the main GO term of which was associated with nucleosome assembly. Considering the function of histones in mitosis, it was concluded that the nucleosome and chromatin assembly may be modified in BC.

MAF encodes a nuclear transcription-regulating protein characterized by a basic region and leucine zipper domain; it has crucial roles in a variety of cellular processes (29). MAFG is a small member of the MAF protein family that consists of little more than the DNA binding and dimerization motif (30), which is able to partially co-localize with FBJ murine osteosarcoma viral oncogene homolog (FOS) in the nucleus and heterodimize with it (31). Members of the AP-1 family of transcription factors are dimeric complexes involved in cellular proliferation, transformation and death (32), and contains the FOS, jun proto-oncogene (JUN) and activating transcription factor (ATF) protein families, of which FOS is a major member. AP-1 family members are immediate early genes induced by a variety of stress signals and control the stress response including cell proliferation, apoptosis and tumorigenesis (33). By forming heterodimers, FOS proteins aid the binding of AP-1 to DNA and exert oncogene activity, leading to tumorigenesis (34). Previous reports demonstrated that the genes encoding AP-1 participate in the cancer-associated immune and inflammation pathways in BC (6). The data from the present study revealed that the expression of MAFG in BC was upregulated, which may increase tumorigenesis by promoting the formation of FOS dimers and the AP-1 complex.

CDC6 is an essential regulator of DNA replication in eukaryotic cells that assembles pre-replicative complexes at origins of replication during the G1 phase of the cell division cycle (35). MCM family members encode highly conserved proteins that act as enzymatically active helicases (36). MCMs drive the formation of pre-replicative complexes (PRCs), the formation of which is the first key event during the G1 phase of the cell cycle (37). MCMs and CDC6 are key proteins in the mechanism of DNA replication and are functionally associated with each other during the cell cycle (38). CDC6 is responsible for the loading of MCM proteins onto origins of replication and, with the presence of CDC6, MCMs bind to the chromatin specifically during the G1 phase of the cell cycle (35). The increased expression of CDC6 and MCM has been observed in dysplastic cells and CDC6 and MCM are consequently considered to be specific biomarkers of proliferating cells (38). Recent studies have revealed the proto-oncogenic activity of CDC6, with its overexpression interfering with the expression of certain tumor suppressor genes and potentially promoting DNA hyper-replication, inducing a senescence response similar to that caused by oncogene activation (39,40). Some members of the MCM family, such as MCM7, were also overexpressed and amplified in a variety of human malignancies (41). In the present study, the expression of CDC6 was found to be upregulated, revealing that in BC cells, DNA replication is aberrant.

In conclusion, the present study reveals that the genes PGR, MAFG, CDC6 and MCMs, and a set of histones are important factors in BC and have key roles in mitotic processes, including nucleosome assembly, spindle checkpoint and DNA replication. However, the small size of the microarray sample and the lack of experimental variation are limitations. Thus, further studies with larger sample size and experimental verification should be performed to confirm the conclusions of the current study.