Ruguo key genes and tumor driving factors identification of bladder cancer based on the RNA-seq profile.

AIM: This study aimed to select several signature genes associated with bladder cancer, thus to investigate the possible mechanism in bladder cancer.
METHODS: The mRNA expression profile data of GSE31614, including ten bladder tissues and ten control samples, was downloaded from the Gene Expression Omnibus. The differentially expressed genes (DEGs) in bladder cancer samples compared with the control samples were screened using the Student's t-test method. Functional analysis for the DEGs was analyzed using the Database for Annotation, Visualization, and Integrated Discovery from the Gene Ontology database, followed by the transcription function annotation of DEGs from Tumor-Associated Gene database. Motifs of genes that had transcription functions in promoter region were analyzed using the Seqpos.
RESULTS: A total of 1,571 upregulated and 1,507 downregulated DEGs in the bladder cancer samples were screened. ELF3 and MYBL2 involved in cell cycle and DNA replication were tumor suppressors. MEG3, APEX1, and EZH2 were related with the cell epigenetic regulation in bladder cancer. Moreover, HOXB9 and EN1 that have their own motif were the transcription factors.
CONCLUSION: Our study has identified several key genes involved in bladder cancer. ELF3 and MYBL2 are tumor suppressers, HOXB9 and EN1 are the main regulators, while MEG3, APEX1, and EZH2 are driving factors for bladder cancer progression.

Introduction

Bladder cancer is one of the most common malignancies in epithelial and mesenchymal tissues of human’s urinary system and contributes to the increasing morbidity and mortality.1 Statistics data refer that approximately 74,000 new cases were diagnosed as bladder cancer in 2014.2 Treatment methods such as surgery, chemotherapy, and drug use play significant role in improving quality of life of patients with bladder cancer.3,4 However, prognosis of bladder cancer is unsatisfactory, with the poor 5-year survival rate due to difficulty in diagnosing and easy metastasis.4 Therefore, it will be of great significance to explore several signature genes and biomarkers for the target treatment and diagnosis of bladder cancer.

Recently, research on the signature genes and transcription factors (TFs), which are closely associated with bladder cancer, is becoming the hot spot. For instance, Dyrskjøt et al5 investigated that 12 gene signatures were associated with nonmuscle invasive bladder cancer prognosis, such as UBE2C, COL18A1, SKAP2, and NEK1. According to Riester et al,6 several tumor metastasis, tumor stage, and progression genes such as FN1, NNMT, POSTN, and SMAD6 are the prognostic signatures for high-risk bladder cancer. Dancik et al7 identified that genetic alterations of TP53, HRAS, KEM6A, and FGFR3 were the diagnostic signature genes that will drive the muscle invasive bladder cancer. Besides, increasing evidences suggest that motifs with the conserved short sequences in the upstream regions of gene promoters play crucial roles in gene transcription process in many cancers, such as bladder cancer,8 prostate cancer,9 and colorectal cancer.10 Rachakonda et al11 reported that mutation of the telomerase reverse transcriptase (TERT) promoter was a potential prognostic target for bladder cancer in clinical application. Despite many studies on the exploration of the prognostic biomarkers in bladder cancer, the mechanism for poor diagnosis and prognoses of bladder cancer still remain largely unknown.

In this study, bioinformatics approaches were used to screen the differentially expressed genes (DEGs) and driving factors that may be helpful for the bladder cancer treatment and prognosis based on the RNA sequencing profile data. The significant functions of the DEGs and the motifs of crucial genes in bladder cancer were analyzed to predict the key genes and driving factors for bladder cancer. This study aimed to explore several signature genes for the target treatment of bladder cancer and to investigate the mechanism of bladder cancer. Our study may provide theoretical basis for the diagnostic and therapeutic exploration for bladder cancer treatment.

Methods

Data resources and preprocessing

The mRNA expression profile data of GSE31614 was downloaded from the Gene Expression Profile in National Center for Biotechnology Information (http://www.ncbi.nlm.nih.gov/geo/) based on the platform of GPL9115 Illumina Genome Analyzer II (Homo sapiens). The study contains the mRNA expression of cancer and matched adjacent tissues of seven testicular germ cell tumors and ten transitional cell carcinomas of bladder. In this study, the ten bladder cancer samples matched with ten matched normal samples were chosen for subsequent analysis.

Reads alignment and transcript expression analysis

The reads of genes sequences in the genome of patients with bladder cancer were aligned with the hg19 genome sequence deposited in the University of California Santa Cruz Genome Browser Database using the Bowtie software.12 Read with only one alignment result and with the base mismatch no more than two was selected. In addition, the expression values of genes in each transcript in each sample were calculated using the Cufflinks tool,13 based on the fragments per kilobase per million mapped fragments method.

Screening of DEGs

The DEGs in bladder cancer samples compared with matched normal samples were identified using the Student’s t-test method14 in R. P<0.05 and |log 2 fold change| ≥2 were chosen as the threshold.

Function annotation of DEGs

The gene ontology functions of the selected DEGs, including biological process (BP) and cellular component, were analyzed using the Database for Annotation, Visualization, and Integrated Discovery software.15 The false discovery rate <0.005 was chosen as the cutoff criterion. Besides, the genes that have the transcription regulatory functions were analyzed and selected based on the information of transcript database. Finally, genes that functions as oncogenes or tumor suppressors were selected based on the information of Tumor-Associated Gene database (http://sourceforge.jp/projects/sfnet_tdt/).

Prediction of the upstream regulatory elements of DEGs

The region in one gene from upstream 1 kb to the downstream 0.5 kb of the transcription start site was considered as the promoter region in this study. The motif of the selected DEGs in the promoter region was selected using the Seqpos.16 P<0.00001 and z score <-7 were chosen as the threshold.

Results

DEGs screening

With P<0.05 and |log 2 fold change| ≥2 as the threshold, a total of 1,571 upregulated and 1,507 downregulated DEGs in the bladder cancer samples were compared with matched normal samples using Student’s t-test method. Besides, 57 upregulated DEGs and 79 downregulated DEGs were analyzed as the TFs based on the transcript database analysis. Three upregulated genes such as E74-like factor 3 (ELF3), v-myb avian myeloblastosis viral oncogene homologue-like 2 (MYBL2), and MYCN were the known oncogenes, while three downregulated genes, including EHF, FNACB, and PML were the known tumor suppressors.

Functional analysis for the DEGs

The selected downregulated DEGs significantly participated in the BP terms such as cell adhesion, extracellular matrix organization, response to endogenous stimulus, and chemical homeostasis (Table 1A). Also, the upregulated genes were mainly involved in BP terms related to cell cycle, such as mitotic cell cycle, cell division, DNA replication, and cell cycle checkpoint (Table 1B).

Table 1

The significantly enriched GO terms by the differentially expressed genes in bladder cancer

Term	Count	P-value	FDR
A: downregulated
GO:0007155 – cell adhesion	128	6.08E-20	1.11E-16
GO:0030198 – extracellular matrix organization	32	4.15E-11	7.55E-08
GO:0001501 – skeletal system development	61	1.62E-10	2.94E-07
GO:0009719 – response to endogenous stimulus	68	3.96E-09	7.21E-06
GO:0009725 – response to hormone stimulus	62	1.67E-08	3.04E-05
GO:0019226 – transmission of nerve impulse	60	1.70E-08	3.10E-05
GO:0007517 – muscle organ development	43	1.71E-08	3.11E-05
GO:0007167 – enzyme linked receptor protein signaling pathway	58	4.49E-08	8.17E-05
GO:0006873 – cellular ion homeostasis	61	8.13E-08	1.48E-04
GO:0035295 – tube development	42	1.71E-07	3.11E-04
GO:0030182 – neuron differentiation	66	4.55E-07	8.28E-04
GO:0042127 – regulation of cell proliferation	101	9.09E-07	1.66E-03
GO:0001944 – vasculature development	44	9.58E-07	1.74E-03
GO:0007242 – intracellular signaling cascade	145	1.62E-06	2.94E-03
GO:0048878 – chemical homeostasis	72	1.76E-06	3.21E-03
B: upregulated
GO:0007049 – cell cycle	123	1.54E-14	2.77E-11
GO:0000278 – mitotic cell cycle	75	2.06E-14	3.74E-11
GO:0000087 – M phase of mitotic cell cycle	54	1.13E-13	2.05E-10
GO:0051301 – cell division	59	3.02E-11	5.47E-08
GO:0006259 – DNA metabolic process	77	1.51E-08	2.73E-05
GO:0006260 – DNA replication	40	1.56E-08	2.82E-05
GO:0000075 – cell cycle checkpoint	23	1.22E-06	2.21E-03
GO:0051726 – regulation of cell cycle	52	1.72E-06	3.11E-03
GO:0007059 – chromosome segregation	21	2.60E-06	4.71E-03

Abbreviations: FDR, false discovery rate; GO, gene ontology; DNA, deoxyribonucleic acid.

In addition, the heat maps of DEGs showed that there were 35 upregulated genes and 138 downregulated genes in tumor tissues that were associated with the cell surface and extracellular matrix (Figure 1). Genes, such as ITGA6, LAMB3, MMP1, TGFA, VEGFA, WNT10A, and WNT7A were all involved in the cancer-related pathways, suggesting these genes might be involved in the bladder cancer progression.

Figure 1

Heat maps of the expression of cell surface and extracellular matrix proteins. Orange side columns stand for the upregulated genes and pink side columns stand for the downregulated genes. Red stripes stand for the bladder cancer tissue, while blue stripes stand for the normal bladder tissue. Yellow signal stands for upregulation, bright blue signal stands for downregulation in bladder cancer tissue, and black signal stands for the expression value between yellow and bright blue signals.

Finally, we selected 12 genes that were associated with the cell epigenetic regulation in bladder cancer tissues (Figure 2). Genes such as maternally expressed 3 (MEG3), PLD6, SMCHD1, and ZFP36 were downregulated, while APEX nuclease multifunctional DNA repair enzyme 1 (APEX1), APOBEC3C, APOBEC3F, ATF7IP, CHEK1, CNOT6, enhancer of zeste homologue 2 Drosophila (EZH2), and GATAD2A were upregulated in the cancer tissues.

Figure 2

Abnormal expression values of genes associated with the apparent regulation in bladder cancer tissues.

Abbreviation: FC, fold change.

Enrichment analysis of motif in TFs

A total of nine genes from the upregulated DEGs and ten genes from the downregulated DEGs were respectively selected, which have the motifs in the upstream sequence based on the motif enrichment analysis of the TFs (Figure 3). We selected the DNA motif for upregulated genes, such as CAATAAAA for homeobox B9 (HOXB9), TAATTA for engrailed homeobox 1 (EN1), and AGGAAGGGAAAGGA for PRDM1 (Figure 3A), while the motif for downregulated genes were TTTCT for NR3C1, TGTTCT for E2F6, and ACATGCGACATG for TP53 (Figure 3B), suggesting their potential tumor driving roles in bladder cancer (Figure 4).

Figure 3

The expression of the genes that has transcription function in bladder cancer.

Notes: (A) Expression of the upregulated genes that has transcription function. (B) Expression of the downregulated genes that has transcription function.

Abbreviation: FC, fold change.

Figure 4

Motif enrichment of transcription factors in bladder cancer.

Discussion

Bladder cancer is one of the most common malignancies of the human urinary system, with increasing morbidity and mortality due to difficulty in diagnosing and easy metastasis. To explore several signature genes and biomarkers for the target treatment, diagnosis of bladder cancer will be of great significance. In this study, we analyzed the mRNA expression profile of bladder cancer to screen the signature genes and key driving factors for bladder cancer diagnosis and target treatment based on the RNA sequencing. The upregulated ELF3 and MYBL2 involved in cell cycle and DNA replication were the tumor suppressors. Downregulated MEG3 and upregulated APEX1 and EZH2 were involved in the cell epigenetic regulation in bladder cancer. Moreover, HOXB9 and EN1 that have their own motifs were the TFs in bladder cancer samples.

ELF3 is a member of the E26 transformation-specific (ETS) family that can bind and transactivate ETS sequences containing the consensus nucleotide core sequence GGA.17 Seth and Watson18 proved that ETS TFs could activate or repress the expression of genes that are involved in the cell proliferation, differentiation, and development and apop-tosis of tumor cells in many cancers. Overexpression of ELF3 induces the endogenous transforming growth factor beta (TGF-β) type II receptor expression in human breast cancer cells,19 and the TGF-β type II receptor is necessary for mediating the effects of TGF-β on tumor cell growth inhibition.20 Thus, ELF3 may be associated with tumor cell growth. Also, ELF3 has been reported as the candidate transcription regulator involved in human urinary cytodifferentiation.21 In this study, the upregulated ELF3 was selected as the tumor suppressor gene, indicating the inhibitory roles in bladder cancer. Meanwhile, MYBL2 is a member of the MYB family, which is a nuclear protein involved in cell cycle progression.22 Sala23 proved that MYBL2 regulated the cell cycle of cancer cells to affect the tumorigenesis. Also, Wu et al24 reported that genetic variation in cell cycle control genes were related to the increased risk of bladder cancer. Our data showed that the upregulated MYBL2 was involved in the cell cycle function in bladder cancer samples, suggesting that MYBL2 might be the tumor suppressor gene for bladder cancer risk through cell cycle.

Our findings showed that the upregulated HOXB9 and EN1 that have their own motifs were the TFs in bladder cancer samples. HOXB9 is a Abd-B homeobox family protein with a homeobox DNA-binding domain on chromosome 17, while EN1 has a role in controlling development during the central nervous system.25,26 However, it has been reported that overex-pression of HOXB9 promoted the cell proliferation and angio-genesis that was related to breast cancer.27 Also, high level of HOXB9 in breast cancer induces the expressions of several angiogenic factors such as interleukin-8, vascular endothelial growth factor, and ErbB to activate their respective pathways, leading to the lung cancer metastasis.28 Besides, Bell et al29 demonstrated that significant hypermethylation of EN1 at the transcriptional start site was observed, suggesting it may be a novel biomarker for poor prognosis of human salivary gland adenoid cystic carcinoma. Therefore, HOXB9 may contribute to tumorigenesis and EN1 may be related to cancer diagnosis. On the other hand, E2F1 is the candidate gene that binds to the promoter region of HOXB9 to induce its overexpression and then accelerated breast cancer progression.30 EN1 has the Pax5 binding site on the promoter region and Pax2 could regulate the enhancer of Pax5 at the midbrain–hindbrain boundary.31,32 Also, inhibited Pax2 results in cell death of prostate cancer.33 Based on our results, we speculate that HOXB9 with the conserved motif CAATAAAA may promote bladder cancer development, while EN1 with DNA sequence TAATTA might be related to bladder cancer diagnosis.

Meanwhile, our results displayed that the downregulated MEG3 and upregulated APEX1 and EZH2 were the genes involved in the cell epigenetic regulation in bladder cancer samples compared with the controls, suggesting their important roles in bladder cancer. Ying et al34 proved that downregulation of MEG3 activated autophagy and increased cell proliferation of bladder cancer by affecting the chromosome. Role of APEX1 and EZH2 in bladder cancer have not been fully reported in previous researches. However, Kim et al35 proved that APEX1 was the driving factor for colon cancer progression, while Varambally et al36 referred that EZH2 was involved in prostate cancer progression. Based on our study, we speculate that MEG3, APEX1, and EZH2 may be the tumor driving factors for bladder cancer progression.

Conclusion

Our findings indicate that ELF3, MYBL2, HOXB9, EN1, MEG3, APEX1, and EZH2 are identified as key genes involved in bladder cancer progression or development. ELF3 and MYBL2 are the key genes that play crucial roles in suppressing bladder cancer, while MEG3, APEX1, and EZH2 are the driving factors for bladder cancer progression. Also, HOXB9 and EN1 are the TFs that play key roles in bladder cancer development and progression. Our study may provide theoretical basis for the future bladder cancer investigation. However, there were still some limitations in the current study. First, the sample size enrolled in our study was small. Second, in vivo and in vitro experimental verification were not performed to validate the functions of crucial genes in bladder cancer samples. Further studies based on a larger sample size and experiments, such as expression validation or knockdown assay in bladder cancer cell lines, are still needed to confirm our results.