Sage Journals: Discover world-class research

Abstract

BACKGROUND:

Accumulating evidence shows that clinical factors alone are not adequate for predicting the survival of patients with urothelial bladder cancer (UBC), and many genes have been found to be associated with UBC prognosis.

PURPOSE:

The objective of this study is to develop a signature which integrates clinical and molecular information to predict the overall survival of UBC patients more accurate.

MATERIALS AND METHODS:

We integrated messenger RNA (mRNA) and microRNA (miRNA) expression profiles and the corresponding clinical data of 402 UBC patients and 19 normal controls from The Cancer Genome Atlas. Univariate Cox regression followed by a multiple testing correction and an elastic net-regulated Cox regression were adopted to identify a prognostic signature.

RESULTS:

We generated an integrated clinical-RNA signature which consisting of 3 clinical variables, 3 protective mRNAs, 7 risky mRNAs, 2 protective miRNAs and 1 risky miRNA. The area under the receiver operating characteristic curve of the integrated clinical-RNA signature was 0.802, larger than that of the clinical-alone signature (0.709) or the RNA-alone signature (0.726). UBC patients in the high-risk group had a significantly shorter overall survival time compared with patients in the low-risk group (clinical-RNA signature, hazard ratio $=$ 2.441).

CONCLUSIONS:

Our conclusions that we have identified an integrated clinical-RNA signature that was superior to the traditional clinical-alone signature for ascertaining the overall survival prognosis of patients with UBC. These findings provide some novel genes for tumor molecular biologist to further study their functions and mechanisms in UBC tumorigenesis and malignance, and may be useful for effective clinical risk management of UBC patients.

Keywords

Urothelial bladder cancer clinical RNA signature prognosis

1. Introduction

Urothelial bladder cancer (UBC) is most common pathological subtype of bladder cancer, and the clinical prognosis of patients with UBC differs significantly which is hard to predict. If a new criterion had been developed to identify the group of high-risk patients ahead, they could get more suitable treatment including extended curative resection or higher radiation/chemo doses, which may greatly ameliorate the poor survival of UBC patients. So, an accurate prediction model is essential for efficient management of UBC which may bring better clinical outcomes to patients with UBC.

A variety of clinical variables in UBC have been studied with regard to their ability to predict disease recurrence, response to therapy, progression, and survival [1, 2, 3, 4, 5]. In the clinical guideline, clinical stage and tumor grade had been recommended as the most important prognostic indicators to forecast disease recurrence and progression [6]. However, significant variability in patient outcome is observed perhaps due to the underlying molecular heterogeneity within the clinically homogeneous tumor groupings.

In recent years, efforts have been made to better understand UBC through basic researches (molecular and genetic), and many studies have revealed that RNAs could serve as novel biomarkers for UBC prognosis [7, 8, 9, 10, 11, 12, 13, 14]. Kim et al. [9] reported that a four-gene signature could predict prognosis of patients with muscle-invasive bladder cancer. Ratert et al. [11] described that expressions of miR-141 and miR-205 were associated with overall survival in patients with bladder cancer. However, these studies mainly focused on messenger RNA (mRNA) profiles or microRNA (miRNA) profiles of UBC genomes independently, and the prognostic signatures were quite different due to a low sample size or the use of an inappropriate regression method for parameter estimation.

Biologically, identification of prognostic signature by integration of mRNA and miRNA expression profiles incorporating more information and may discover more robust prognostic signature due to potential interactions between miRNAs and mRNAs. Methodologically, in genomic expression analysis, there is a so-called “curse of dimensionality” problem in that the number of genes is much larger than the number of samples [15]. In this setting, ordinary regression is subject to over-fitting and instable coefficients, and stepwise variable selection methods do not scale well [16, 17]. Regression by penalization methods has been successfully adapted to high-dimensional multiple genomic datasets and outperforms univariate and multivariate regression methods [18]. At present, the most commonly used penalization methods are ridge regression, Least Absolute Shrinkage and Selection Operator (LASSO) regression and a hybrid of these (elastic net regression); all three methods are based on penalizing the L1 norm, the L2 norm, and both the L1 norm and L2 norm with tuning parameters. Although the traditional Cox proportional hazards model is widely used to discover cancer prognostic factors, it is not appropriate for the genomic setting due to the high dimensionality and colinearity. Several groups have proposed to combine the Cox regression model with the elastic net dimension reduction method to select survival-correlated genes within a high-dimensional expression dataset and have made available the associated computation procedures [19, 20, 21].

By characterizing genetic alterations, epigenetic alterations, and the expression of cancer genomes, The Cancer Genome Atlas (TCGA) project has provided a comprehensive way to understand UBC [22]. As the number of UBC samples in this project grows, the opportunity to identify prognostic molecular signature for patients with UBC is increasing. In the current study, we subjected the integrated mRNA and miRNA profiles of UBC patients to elastic net-regulated Cox regression analysis and identified an integrated prognostic signature that can predict the overall survival (OS) of UBC patients. These may be helpful for selecting high risk UBC patients for better clinical management.

2. Materials and methods

2.1 Patient characteristics and integrated RNA profiles

Clinical data and RNA expression profiles of UBCs were retrieved from TCGA (https://cancergenome.nih.gov/). The dataset contained 412 UBC patients and 19 patients have matched normal samples. A total of 402 patients were included for further study, and 10 patients were excluded from the study by the following criterions: First, the patients that not having both fully characterized mRNA and miRNA profiles. Second, the patients were alive but the last followup times were unavailable. mRNA expression data were generated by RNA-Seq with HTSeq-Counts workflow type and miRNA expression data were generated by miRNA-Seq with GCGSC miRNA profiling workflow type. The level-3 raw counts of mRNA and miRNA sequencing data were downloaded, combined, full-quantile-normalized, and quantile-filtered.

2.2 Prognostic clinical factors

Survival analysis with log-rank test was preliminarily applied to filter out the demographic and clinical factors that were correlated with the OS of UBC patients from among 9 factors: gender, age at initial diagnosis, clinical stage, tobacco use, anatomic neoplasm subdivision, Karnofsky performance status, human race, lymphovascular invasion, and neoplasm histologic grade. Demographic and clinical factors resulting in a log-rank test P value $<$ 0.05 were further selected into a multivariate stepwise Cox regression analysis to evaluate the contribution of the preliminary selected clinical factors as independent clinical prognostic factors.

2.3 Prognostic mRNAs and miRNAs

mRNAs and miRNAs that differentially expressed between UBCs and normal samples were selected using “glmLRT” method implemented in “edger” R package with FDR $\leqslant$ 0.01 and logFC $\geqslant$ 1 as parameters. A univariate Cox regression method was used to assess the relationship between OS and expressions of differentially expressed RNAs. To alleviate the effect of “curse of dimensionality” for general multivariate Cox regression analysis, we adopted an elastic net-regulated Cox regression analysis with 200,000 iterations and 10 cross-validations [20, 21] to estimate the regression coefficients more accurate. RNAs with elastic net-regulated Cox regression coefficients $\neq$ 0, which were considered candidate RNAs, were included in the RNA signature. RNAs that had a hazard ratio (HR) for death $>$ 1 were considered to be risky RNAs, and those with an HR $<$ 1were defined as protective RNAs.

2.4 Prognosis index and integrated clinical-RNA signature

We adopted a prognosis index (PI) which defined as a linear combination of the observed value weighted by the regression coefficients as indicator for UBC patient prognosis. Specifically,

${\mathrm{PI}}_{j}=\sum\nolimits_{i}{\beta_{i}\times G_{ij}}$

where, PI ${}_{j}$ is the prognostic index of the jth observation, $\beta_{i}$ is the regression coefficient of the $i$ th variable. $G_{ij}$ is the observed value of the $i$ th variable in the $j$ th observation. For the RNA signature, $G_{ii}$ is the relative expression value of the $i$ th RNA in the $j$ th patient, $\beta_{i}$ is the elastic net-regulated Cox regression coefficient of the $i$ th RNA, and the calculated PI was termed as RNA-PI. For the clinical factors, $G_{ij}$ is the observed value of the $i$ th clinical factor in the $j$ th patient, $\beta_{i}$ is the general multivariate Cox regression coefficient of the $i$ th clinical factor, and the calculated PI was termed as clinical-PI.

We also used multivariate stepwise Cox regression analysis to determine the prognostic contributions and independence of RNA-PI and the associated clinical factors [15, 23]. RNA-PI and clinical factors with a multivariate Cox regression P value $<$ 0.05 were considered to be independent predictors for UBC patient prognosis. Furthermore, we constructed a new PI termed as clinical-RNA-PI as an integrated indicator of RNA expressions and clinical factors for UBC patient prognosis. The clinical-RNA-PI was also calculated as a linear combination. For both forms of PI, a weighed prognostic index (WPI) was calculated as the standard form of the PI. i.e.,

$\mathrm{WPI}=\frac{\textit{PI}-\textit{mean}(\textit{PI})}{sd(\textit{PI})}$

where, Mean(PI) and sd(PI) are the mean and standard deviation of the PI vector, respectively.

2.5 Model evaluation

UBC patients that have WPI $>$ 0 were classified into a high-risk group and the others were classified into a low-risk group. OS between the high-risk group and the low-risk group was compared using Kaplan-Meier method with 2-sided log-rank test. The ability and efficiency of the integrated clinical-RNA signature to predict UBC patient outcome was assessed by calculating the area under the curve (AUC) of the time dependent receiver operating characteristic (ROC). Analysis was conducted by the “survivalROC” package in R software for time dependent ROC analysis at 5 years of OS time [24]. Furthermore, we compared the prognostic ability that indicated by the ROC curves of the integrated clinical-RNA signature with the RNA-alone signature and the clinical-alone signature.

2.6 GO and pathway of mRNA and targets of miRNAs

Validated mRNA and targets of miRNAs were available from miRTarBase [25]. Only miRNA and mRNA target interactions with strong evidences, either by reporter assay or western blot, were extracted. Gene ontology enrichment and pathway enrichment were carried out based on these validated targets using online tool STRING [26, 27, 28]. The interactions among products of these targets were visualized by the same tool.

Table 1
Univariate survival analysis and multivariate Cox regression of demographic and clinical variables

Variables	Deaths/Patients (%)	Median	95% CI	Log-rank	Multivariate Cox	Multivariate
		survival		test P	HR (95%) CI	Cox $P$
Age at diagnosis
$\leqslant 65$ y	45/149 (30.20)	2828	1423–NA	9.24E $-$ 05	2.26 (1.44–3.56)	4.34E $-$ 04
$>$ 65 y	131/252 (51.98)	819	602–1064
Clinical stage
Stage I	0/2 (0)	–	–	3.71E $-$ 07	1.47 (1.11–1.95)	7.43E $-$ 3
Stage II	35/128 (27.34)	2828	1869–NA
Stage III	56/138 (40.58)	1367	795–NA
Stage IV	84/132 (63.64)	593	508–739
Lymphovascular invasion
No	39/114 (34.21)	2828	1556-NA	2.12E $-$ 05	1.83 (1.18–2.86)	7.34E $-$ 3
Yes	69/125 (55.20)	623	547–1348
Gender
Female	51/104 (49.04)	941	623–NA	0.364
Male	125/298 (41.95)	1077	819–1804
Tobacco use
Non-smoker	5/13 (38.46)	778	254–NA	0.573
Current smoker	88/220 (40.00)	1163	823–2828
Current reformed smoker	83/169 (49.11)	974	685–1869
Anatomic neoplasm subdivision
Bladder-NOS	89/204 (43.63)	1367	1004–2641	0.379
Dome	11/24 (45.83)	712	547–NA
Neck	1/7 (14.29)	–	–
Trigone	13/31 (41.94)	795	623–NA
Wall anterior	8/24 (33.33)	835	617–NA
Wall lateral	19/46 (41.30)	904	577–NA
Wall NOS	12/22 (54.55)	510	415–NA
Wall posterior	19/37 (55.89)	706	579–NA
Karnofsky performance status
$\leqslant 60$	3/14 (21.43)	630	630–NA	0.559
$\leqslant 80$	15/30 (50.00)	734	612–NA
$\leqslant 100$	43/81 (53.01)	1670	823–2828
Neoplasm histologic grade
High grade	174/378 (46.03)	1005	778–1670	0.119
Low grade	2/21 (9.52)	904	–
Human race
Asian	9/44 (20.45)	1670	674–NA	0.262
Black ${}^{*}$	13/23 (56.52)	696	510–NA
White	147/318 (46.23)	1005	778–1718

Abbreviation: CI, confidence interval. P-values were obtained from log-like test for survival analysis and Wald test for Cox regression, patients were omitted when data is unavailable. For multivariate Cox regression, variable coding is as follows: Age at diagnosis (1 $\leqslant$ 65 years; 2 $\geqslant$ 65 years), Clinical stage (1 $=$ leve I; 2 $=$ level II; 3 $=$ level III; 4 $=$ level IV), Lymphovascular invasion (0 $=$ No; 1 $=$ Yes).

2.7 Statistical analysis

All statistical analyses were completed by R software [29] (version 3.3.2). Survival analysis with log-rank test was used to compare survival distributions between different groups, and it is implemented by “survfit” function that incorporated in “survival” R package. Univariate Cox regression analysis was used to evaluate the association between RNA expressions across samples with OS of UBC patients. Multivariate Cox regression analysis with step-wise procedure to select independent variables. Both univariate and multivariate Cox regression analyses with Wald test were implemented by “coxph” function that incorporated in “survival” R package. The generalized linear model (GLM) likelihood ratio test is based on the idea of fitting negative binomial GLMs with the Cox-Reid dispersion estimates. The testing can be done by using the “glmLRT” function implemented in “edger” R package [30]. The elastic-net penalized Cox regression was implemented by “glmnet” function and cross validation was implemented by “cv.glmnet” function incorporated in “glmnet” R package. All statistical tests were considered to be significant with $P<$ 0.05 or FDR $<$ 0.05, as appropiate.

Figure 1.

Cross-validation error curve. The left vertical line shows where the cross-validation error curve hits its minimum (lambda $=$ 0.05). The right vertical line shows the most regularized model with cross-validation error within 1 standard deviation of the minimum. The minimum was achieved by a fairly regularized model ( $n=$ 13), but the right line indicates that the null model (no coefficients included). The numbers at the top of the figure indicate the number of the nonzero coefficients.

Table 2

RNAs correlated with UBC prognosis

No	miRNA	P Value ${}^{\$}$	FDR	Type ${}^{\#}$
1	hsa-mir-200c	1.11E-04	2.32E-02	Protective
2	hsa-mir-598	1.35E-04	2.32E-02	Protective
3	hsa-mir-143	3.24E-04	3.71E-02	Risky
4	RUNX2	7.09E-05	8.61E-04	Risky
5	LAMA2	1.67E-04	1.37E-02	Risky
6	UBD	2.20E-04	1.52E-02	Protective
7	INHBB	2.59E-04	1.65E-02	Risky
8	CATSPER2	5.42E-04	2.18E-02	Protective
9	KCNK6	5.81E-04	2.26E-02	Risky
10	PLSCR4	6.02E-04	2.27E-02	Risky
11	ZNF600	1.83E-03	3.94E-02	Protective
12	PAM	2.08E-03	4.05E-02	Risky
13	PDGFD	2.21E-03	4.16E-02	Risky

${}^{\$}$ P-values were obtained from Wald test in univariate Cox’s model. ${}^{\#}$ Types included protective (HR $<$ 1) and risky (HR $>$ 1).

3. Results

3.1 Clinical characteristics correlated with UBC patient survival

Table 3
Multivariate Cox regression of demographic and clinical variables

Variables	HR	95% CI	Multivariate Cox P
RNA-PI	2.351	1.565–3.533	3.87E-05
Age at diagnosis	2.214	1.404–3.491	6.29E-04
Clinical stage	1.382	1.053–1.813	0.0195
Lymphovascular	1.690	1.102–2.593	0.0163
invasion
			1.907E-12

P-values of multivariate Cox regression model were obtained from Wald test.

Of the 402 patients in the TCGA UBC cohort, 176 patients were deceased and 226 were alive at the time of last follow-up. The average age at initial diagnosis of the patients in this cohort was 60.14 years, and the median overall survival time was 813 days (95% confidence interval [CI], 731–895 days). Demographic and clinical data for the TCGA UBC cohort were summarized in Table 1. Gender, age at initial diagnosis, clinical stage, tobacco use, anatomic neoplasm subdivision, Karnofsky performance status, human race, lymphovascular invasion, and neoplasm histologic grade were evaluated by survival analysis. Among these clinical characteristics, only age at diagnosis, lymphovascular invasion, and clinical stage had a P value less than 0.05 by Log-rank test (Table 1). Multivariate Cox regression analysis of these 3 factors revealed that all of them were independent prognostic factors (Table 1). The results of this preliminary assessment demonstrated that the survival data in the TCGA UBC cohort, despite the censored data, were informative and applicable to identifying a prognostic signature.

Figure 2.

WPI distributions and stratifications of UBC patients of clinical signature (A), RNA signature (B) and clinical-RNA signature (C), respectively.

Figure 3.

Kaplan-Meier plots for stratifications based on clinical signature (A), RNA signature (C) and clinical-RNA signature (E), respectively. ROC for WPIs based on clinical signature (B), RNA signature (D) and clinical-RNA signature (F), respectively.

Table 4

Slim gene ontology enrichment

Category	ID	Pathway description	Observed gene count	FDR
Biological	GO.0016925	Protein sumoylation	13	5.9E $-$ 19
processes	GO.0018205	Peptidyl-lysine modification	13	6.12E $-$ 13
	GO.0043687	Post-translational protein modification	13	5.45E $-$ 12
	GO.0070647	Protein modification by small protein conjugation or removal	14	6.59E $-$ 09
	GO.0018193	Peptidyl-amino acid modification	14	1.91E $-$ 08
	GO.0045892	Negative regulation of transcription, DNA-templated	15	3.42E $-$ 08
	GO.0010629	Negative regulation of gene expression	16	5.68E $-$ 08
	GO.0032446	Protein modification by small protein conjugation	12	1.75E $-$ 07
	GO.1903507	Negative regulation of nucleic acid-templated transcription	14	4.01E $-$ 07
Cellular	GO.0031519	PcG protein complex	7	3.58E $-$ 11
components	GO.0035102	PRC1 complex	5	2.99E $-$ 09
	GO.0016604	Nuclear body	10	2.97E $-$ 08
	GO.0016605	PML body	6	1.69E $-$ 06
	GO.0043234	Protein complex	21	1.69E $-$ 06
	GO.0044451	Nucleoplasm part	10	8.03E $-$ 06
	GO.0000803	Sex chromosome	4	1.19E $-$ 05
	GO.0001739	Sex chromatin	3	1.39E $-$ 05
	GO.0044427	Chromosomal part	9	6.79E $-$ 05
Molecular	GO.0019789	SUMO transferase activity	5	8.20E $-$ 09
function	GO.0005515	Protein binding	24	4.97E $-$ 07
	GO.0003682	Chromatin binding	10	2.45E $-$ 06
	GO.0071535	RING-like zinc finger domain binding	3	3.93E $-$ 06
	GO.0008134	Transcription factor binding	8	1.15E $-$ 04
	GO.0003713	Transcription coactivator activity	7	2.05E $-$ 04
	GO.0044877	Macromolecular complex binding	11	2.05E $-$ 04
	GO.1901363	Heterocyclic compound binding	22	2.84E $-$ 04
	GO.0097159	Organic cyclic compound binding	22	3.22E $-$ 04

3.2 Prognostic RNA signature for UBC patients

Differential expression analysis revealed that 3043 RNAs were differentially expressed between UBCs and normal samples. Among these RNAs, 211 RNAs (adjusted P value $<$ 0.05) were preliminarily identified as OS predictors by univariate Cox regression. Finally, elastic net regulated Cox regression resulted in the selection of 10 mRNAs and 3 miRNA (Table 2) as the prognostic RNA signature with minimal lambda $=$ 0.05 (Fig. 1). Multivariate stepwise Cox regression analysis revealed that RNA-PI, age at initial diagnosis, clinical stage, and lymphovascular invasion were independent prognostic predictors for UBC patient survival (Table 3). The HR (HR $=$ 2.351) of the RNA signature was greater than that of the demographic and clinical variables (Table 3), and it implied that the RNA signature had superior performance compared with traditional clinical variables. We also calculated PIs for clinical-alone and clinical-RNA signatures. As continuous variables, the RNA-PI and clinical-RNA-PI were both significantly correlated with UBC patient survival by univariate Cox regression analysis (HR $=$ 18.12; 95% CI $=$ 8.361–39.29; $P=$ 2.136E-13 for RNA-PI; HR $=$ 3.113; 95% CI $=$ 2.304–4.206; $P=$ 1.37E-13 for clinical-RNA-PI).

The WPIs ranged from $-$ 2.13 to 1.46 for clinical-signature (Fig. 2A), $-$ 5.77 to 3.99 for RNA-signature (Fig. 2B) and $-$ 3.14 to 2.09 for clinical-RNA-signature (Fig. 2C). We split the UBC patients into high-risk and low-risk groups by WPIs for patients that have full age at diagnosis, lymphovascular invasion, and clinical stage information if necessary ( $n=$ 238 for clinical-signature and clinical-RNA signature, $n=$ 402 for RNA-signature). Specifically, patients with WPIs $>$ 0 were classified as high-risk group and patients with WPIs $\leqslant$ 0 were classified as low-risk group (Fig. 2). Survival analysis showed that the median OS time was significantly shorter in the high-risk group than in the low-risk group for all WPI based stratifications (761.2 versus 1007.2 days, HR $=$ 1.794, 95% CI $=$ 1.457–2.208, $P<$ 0.001, for clinical-signature, Fig. 3A; 748.4 versus 905.8 days, HR $=$ 1.876, 95% CI $=$ 1.586–2.220, $P<$ 0.01, for RNA-signature, Fig. 3C; 765.7 versus 1039.2 days; HR $=$ 2.441; 95% CI $=$ 1.927–3.092; $P<$ 0.001, for clinical-RNA-signature, Fig. 3E).

In all of the survival analyses, fewer events occurred after 5 years, we tested the ability of the integrated RNA signature to predict the survival outcome of UBC patients at, and around, this time point by time dependent ROC analysis. The AUCs were 0.709, 0.726 and 0.802 at 5 years of OS for clinical-signature (Fig. 3B), RNA-signature (Fig. 3D) and clinical-RNA signature (Fig. 3F), respectively. Thus, the integrated clinical-RNA signature had superior performance compared with clinical alone and RNA alone signatures whatever in terms to greater HR or greater AUC.

Figure 4.

Network of proteins. A. The network predicted with methods of neighborhood, gene fusion, co-occurrence, co-expression, experiments, databases and text-mining. This network was enriched in interaction with a p-value of 4.22e-15. B. The network predicted only with experiments. Disconnected nodes were hidden. This network was enriched in interaction with a p-value of 0.0324.

Table 5

KEGG pathway enrichment

ID	Term	Observed	FDR
		gene count
5206	MicroRNAs in cancer	6	2.86E-05
3013	RNA transport	5	5.91E-04
4722	Neurotrophin signaling pathway	4	3.27E-03
4120	Ubiquitin mediated proteolysis	4	4.46E-03
4010	MAPK signaling pathway	4	3.93E-02

3.3 Gene ontology and pathway enrichment

STRING database (http://string-db.org/) [28] is employed to analyze the pathway enrichment of miRNAs’ targets (Supplementary Table 1) and mRNAs of the RNA signature. The gene ontology included here are just slim gene ontology, which gives an overview without the details of the specific terms [26]. The results show that these genes are mainly enriched in protein sumoylation biological process (Table 4). Others are enriched in transcription biological process. Protein sumoylation is closely associated with interaction and regulation to protein-protein interaction in living cells. And it is also closely associated with tumorigenesis [31].

In cellular component analysis, we find that many gene products locate in chromosome and involve in regulating target genes by chromosome modification (Table 4). These gene products generally constitute some protein complex to bind chromatin.

In molecular function analysis (Table 4), molecular function of SUMO transferase activity is a major function in UBC patients. The enrichment of molecular function supports biological process. Besides, we can find that many genes are involved in protein binding, chromatin binding and domain binding. The main functions of these genes were binding and activity. The results indicated that these genes as direct regulators for regulation gene expression.

In KEGG pathway analysis, miRNA’s target genes and mRNA are mainly enriched in four KEGG pathways (FDR $<$ 0.01) (Table 5). The enrichment of KEGG has some different from biological process. They are miRNAs in cancer, RNA transport, neurotrophin signaling pathway, and ubiquitin mediated proteolysis. In these four pathways, RNA transport, ubiquitin mediated proteolysis, and MAPK signaling pathway are related with tumor. Of these KEGG pathways, MAPK play critical roles in tumor formation and metastasis.

3.4 Network analysis

STRING was used to visualize the protein-protein interaction network among miRNAs’ targets and mRNAs of the RNA signature. Figure 4A showed the network with methods of neighborhood, gene fusion, co-occurrence, co-expression, experiments, databases and text-mining. This network wags enriched in interaction, which assumed that these protein worked altogether and joined in the related pathways. After removing evidences with lower confidence and only using data experiments, a network with fewer genes was built (Fig. 4B). In both networks, TP53 was the center and had high degrees, indicating the important role of it in the prognosis of UBC.

4. Discussion

In this study, we proposed a novel integrated mRNA, miRNA and clinical signature that could predict UBC patient survival more accurate than clinical-only and RNAs-only signatures. Pathway analysis revealed that genes in our signature were involved in cell death and survival. All results suggested that the integrated signature was suitable and had superior prognostic value for UBC patient prognosis. Our study provides novel insights into the significance of bringing molecular and clinical information together for UBC prognosis.

Although markers for classifying UBC molecular subtypes have been identified, markers associated mainly with the pathogenesis of UBC may not predict the prognosis [32]. Furthermore, it is impertinent to use the clinical characteristics alone to distinguish the high-risk patients due to the underlying molecular heterogeneity within the clinically homogeneous tumor groupings. Thus we objected to identify common RNAs that consistently drive the outcome of UBC patients irrespective of the clinical or molecular subtypes. For UBC prognostic markers discovery, studies have been focused on either mRNA [9, 33] or miRNA profiles [7, 10, 11, 12, 13, 14] with limited sample size. To our knowledge, there is no study have been proposed a combination signature of mRNA expression, miRNA expression, and clinical characteristics to improve the UBC patient prognosis.

Of the prognostic RNAs, hsa-mir-200c was previously reported to be associated with UBC patient survival [34]. It has been reported that hsa-mir-200c regulates epithelial to mesenchymal transition (EMT) and restores expression of E-cadherin in breast and ovarian cancer cells [35, 36, 37]. Yu et al. reported that high levels of hsa-mir-200c expression inhibit cancer invasion and stimulate cancer cell proliferation, possibly via up-regulation of E-cadherin, and that high levels of hsa-mir-200c expression correlate with better survival of patients with curative resection of pancreatic cancer [38]. Furthermore, to the best of our knowledge, we are the first to report the hsa-mir-598, UBD, CATSPER2 and ZNF600 were also associated with UBC patient survival. Expression of miR-598-3p was down-regulated in the serum of breast cancer patients, using miR-598-3p as a biomarker, the sensitivity and specificity for the detection of breast cancer was 95.0 and 87.5% [39]. High expression of UBD correlates with epirubicin resistance and indicates poor prognosis in triple-negative breast cancer [40]. ZNF600 were down regulated in smokers compared to non-smokers [41], and cigarette smoking was an important risk factor for bladder cancer in both sexes [42]. Although the specific function of ZNF600 is unclear, ZNF600 expression is closely associated with smoking. Consistent with the GO result, UBD is not only an ubiquitin protein associated with ubiquitylation, which plays a vital role in survival [43]. CATSPER2 plays a central role in calcium-dependent physiological responses essential for successful fertilization, such as sperm hyperactivation, acrosome reaction and chemotaxis towards the oocyte. The relationship between prognosis of UBC and CATSPER2 expression has not been reported yet. ZNF600 belongs to the krueppel C2H2-type zinc-finger protein family which involved in transcriptional regulation. Transcriptional regulation factor such as ZNF600 generally affects cancer by promoting cell proliferation.

Among the risky genes, hsa-mir-143, RUNX2,LAMA2, INHBB, KCNK6, PLSCR4, PAM andPDGFD are associated with tumor development. In previous report of study, microRNA-143 as a tumor suppressor for bladder cancer and microRNA-143 inhibits cell migration and invasion by targeting matrix metalloproteinase 13 in prostate cancer. Thus, the function of microRNA-143 is heterogeneous in different study. And this condition needs to be further study. RUNX2 and TP53 might be functionally related and are likely involved in bladder tumor carcinogenesis and aggressiveness, and RUNX2 and p53 independently predict early tumor recurrence in bladder carcinoma patients, with the highest prediction accuracy being achieved on their combined high expression [44]. TP53-induced miR-34a contributes liver regeneration termination via regulation of INHBB and MET. Other risky genes have not been reported in bladder cancer yet. However, LAMA2, INHBB and PDGFD all involved in cell migration and embryonic development [45, 46]. In fact, cell migration is closely associated with cancer. On the other hand, embryonic development and cancer development have some similar biological processing such as rapid cell proliferation and EMT pathway. KCNK6, PLSCR4 and PAM performed their special role in bladder cancer. Although these three genes are high-risk genes in bladder cancer, it is difficult to elaborate their functions. According to main functions of PLSCR4 and PAM, they are ion binding proteins. And KCNK6 is a transport for K ${}^{+}$ . Therefore, we inferred that the three genes might be associated with drug resistance. And we expected that these genes can be validated by further biological experiment.

In conclusion, we have identified a 3-clinical and 13-RNA integrated signature that can predict the survival outcome of UBC patients more accurate than RNA-alone or clinical-alone signatures. Our findings may help researchers understanding of UBC cell death and survival, develop targeted therapy, and identify high-risk UBC patients for better disease management.

Footnotes

Acknowledgments

which made the genomic data of UBC available. We also thank the anonymous reviewers for helpful comments. This work was supported by The Postdoctoral Science Foundation of Central South University.

Conflict of interest

The authors have declared that no competing interests exist.

Supplement material

Supplementary Table 1

Mature sequences of miRNA signature and experimentally validated mRNA targets

Mature miRNA	Targets	Mature miRNA	Targets
hsa-mir-143-3p	KRAS	hsa-miR-200c-3p	TUBB3
	MYO6		BMI1
	DNMT3A		SIP1
	FNDC3B		BAP1
	MAPK7		ZEB2
	FSCN1		ZEB1
	HK2		Zeb2
	SERPINE1		Zeb1
	MACC1		FN1
	PTGS2		ZFPM2
	JAG1		UBE2I
	AKT1		PTPN13
	COL1A1		RNF2
	HRAS		RCOR3
	MDM2		BRD7
	BCL2		ACVR2B
	MMP13		MSN
	SDC1		NTRK2
	RREB1		ERRFI1
	CD44		CCNE2
	KLF5		XIAP
	BRAF		BCL2
	TNF		TIMP2
	NR2C2		FBLN5

Supplementary Table 1, continued
Mature miRNA	Targets	Mature miRNA	Targets
	IL13RA1		VEGFA
	DYT10		NCAM1
	COL3A1		Vldlr
	DDX6		Reln
	LIMK1		IKBKB
	HNF4A		FLT1
			KLF9
hsa-miR-143-5p	COX2	hsa-miR-200c-3p	TBK1
	IGF1R		PMAIP1
			NTF3
			LPAR1
			EDNRA
			RHOA
			KLHL20
			PTPRD
			ELMO2
			ERBB2IP
			WDR37
			VAC14
			TCF7L1
			RASSF2
			HOXB5
			RIN2
			KLF11
			SEPT7
			SHC1
			MYB
			ETS1
			USP25
			EFNA1
			RND3
			DNMT3A
			DNMT3B
			SP1
			CFL2
			CDH11
			SEC23A
			KDR
			HFE
			DLC1
			ATRX
			ZNF217
			BTC
			ZFPM1
			PIN1
			KRAS
			NOTCH1
			GATA4
			SUZ12
			ROCK2
			UBQLN1
			E2F3

References

Dalesio

Schulman

C.C.

Sylvester

et al., Prognostic factors in superficial bladder A study of the european organization for research on treatment of cancer: Genitourinary tract cancer cooperative group, Journal of Urology 129(4) (1983), 730–733.

Fitzpatrick

J.M.

West

A.B.

Butler

M.R.

et al., Superficial bladder tumors (stage pTa, grades 1 and 2): The importance of recurrence pattern following initial resection, Journal of Urology 135(5) (1986), 920–922.

Parmar

M.K.B.

Freedman

L.S.

Hargreave

T.B.

et al., Prognostic factors for recurrence and followup policies in the treatment of superficial bladder cancer: Report from the british medical research council subgroup on superficial bladder cancer (urological cancer working party), Journal of Urology 142(2) (1989), 284–288.

Soloway

M.S.

, The relationship among multiple recurrences, progression and prognosis of patients with stages TA and T1 transitional cell cancer of the bladder followed for at least 20 years – editorial comment, Journal of Urology 153(6) (1995), 1826–1827.

Sylvester

R.J.

van der Meijden

A.P.M.

Oosterlinck

et al., Predicting recurrence and progression in individual patients with stage Ta T1 bladder cancer using EORTC risk tables: A combined analysis of 2596 patients from seven EORTC trials, European Urology 49(3) (2006), 466–477.

Hall

M.C.

Chang

S.S.

Dalbagni

et al., Guideline for the management of nonmuscle invasive bladder cancer (stages Ta, T1, and Tis): 2007 update, Journal of Urology 178(6) (2007), 2314–2330.

Rosenberg

Baniel

Spector

et al., Predicting progression of bladder urothelial carcinoma using microRNA expression, Bju International 112(7) (2013), 1027–1034.

Konety

B.R.

, Molecular markers in bladder cancer: A critical appraisal, Urologic Oncology 24(4) (2006), 326–337.

Kim

W.J.

Kim

S.K.

Jeong

et al., A four-gene signature predicts disease progression in muscle invasive bladder cancer, Molecular Medicine 17(5–6) (2011), 478–485.

10.

Y.Q.

Y.Z.

et al., MicroRNA expression profiles in muscle-invasive bladder cancer: Identification of a four-microRNA signature associated with patient survival, Tumor Biology 36(10) (2015), 8159–8166.

11.

Ratert

Meyer

H.A.

Jung

et al., miRNA profiling identifies candidate miRNAs for bladder cancer diagnosis and clinical outcome, Journal of Molecular Diagnostics Jmd 15(5) (2013), 695–705.

12.

Wang

Zhang

Wang

et al., MicroRNA-214 suppresses oncogenesis and exerts impact on prognosis by targeting PDRG1 in bladder cancer, Plos One 10(2) (2015), e0118086.

13.

Fujii

Tatsumi

Fujimoto

et al., Microrna-145 promotes differentiation in human urothelial carcinoma through down-regulation of syndecan-1, BMC Cancer 15(1) (2015), 818.

14.

Zhou

Tang

Xiao

et al., A panel of eight-miRNA signature as a potential biomarker for predicting survival in bladder cancer, Journal of Experimental & Clinical Cancer Research 34(1) (2015), 1–12.

15.

Yoshihara

Tajima

Yahata

et al., Gene expression profile for predicting survival in advanced-stage serous ovarian cancer across two independent datasets, Plos One 5(3) (2010), 1191–1204.

16.

H.F.

Lee

K.L.

and Mark

D.B.

, Multivariable prognostic models: Issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors, Statistics in Medicine 15(4) (1996), 361–387.

17.

Yuan

and Lin

, Model selection and estimation in regression with grouped variables, Journal of the Royal Statistical Society 68(1) (2006), 49–67.

18.

Bøvelstad

H.M.

Nygård

Størvold

H.L.

et al., Predicting survival from microarray data – a comparative study, Bioinformatics 23(16) (2004), 2080–2087.

19.

Simon

Friedman

Hastie

et al., Regularization paths for Cox’s proportional hazards model via coordinate descent, Journal of Statistical Software 39(5) (2011), 1–13.

20.

Tibshirani

Bien

Friedman

et al., Strong rules for discarding predictors in lasso-type problems, Journal of the Royal Statistical Society 74(2) (2010), 245–266.

21.

Friedman

Hastie

and Tibshirani

, Regularization paths for generalized linear models via coordinate descent, Journal of Statistical Software 33(1) (2010), 1.

22.

Network TCGA, Comprehensive molecular characterization of urothelial bladder carcinoma, Nature 507(7492) (2014), 315–322.

23.

Xiong

Bing

et al., An integrated mRNA and microRNA expression signature for glioblastoma multiforme prognosis, Plos One 9(5) (2014), e98419.

24.

Heagerty

P.J.

Lumley

and Pepe

M.S.

, Time-dependent ROC curves for censored survival data and a diagnostic marker, Biometrics 56(2) (2000), 337–344.

25.

Hsu

S.D.

Tseng

Y.T.

Shrestha

et al., miRTarBase update 2014: An information resource for experimentally validated miRNA-target interactions, Nucleic Acids Research 42(1) (2014), 78–85.

26.

Consortium TGO, Gene ontology consortium: Going forward, Nucleic Acids Research 43(1) (2015), 1049–1056.

27.

Kanehisa

Goto

Sato

et al., Data, information, knowledge and principle: Back to metabolism in KEGG, Nucleic Acids Research 42(1) (2014), 199–205.

28.

Szklarczyk

Franceschini

Wyder

et al., STRING v10: Protein-protein interaction networks, integrated over the tree of life, Nucleic Acids Res 43(1) (2015), 447–452.

29.

Ihaka

and Gentleman

, R: A language for data analysis and graphics, Journal of Computational & Graphical Statistics 5(3) (1996), 299–314.

30.

Robinson

M.D.

McCarthy

D.J.

and Smyth

G.K.

, edgeR: A bioconductor package for differential expression analysis of digital gene expression data, Bioinformatics 26(1) (2010), 139–140.

31.

Masoumi

K.C.

Marfany

et al., Putative role of SUMOylation in controlling the activity of deubiquitinating enzymes in cancer, Future Oncology 12(4) (2016), 565–574.

32.

Hurst

C.D.

Platt

F.M.

Taylor

C.F.

et al., Novel tumor subgroups of urothelial carcinoma of the bladder defined by integrated genomic analysis, Clinical Cancer Research 18(21) (2012), 5865–5877.

33.

Riester

Taylor

J.M.

Feifer

et al., Combination of a novel gene expression signature with a clinical nomogram improves the prediction of survival in high-risk bladder cancer, Clinical Cancer Research an Official Journal of the American Association for Cancer Research 18(5) (2012), 1323–1333.

34.

Liu

Qiu

Tan

et al., miR-200c inhibits invasion, migration and proliferation of bladder cancer cells through down-regulation of BMI-1 and E2F3, Journal of Translational Medicine 12(1) (2014), 305.

35.

Hurteau

G.J.

Carlson

J.A.

Roos

et al., Stable expression of miR-200c alone is sufficient to regulate TCF8 (ZEB1) and restore E-cadherin expression, Cell Cycle 8(13) (2009), 2064–2069.

36.

Park

S.M.

Gaur

A.B.

Lengyel

et al., The miR-200 family determines the epithelial phenotype of cancer cells by targeting the E-cadherin repressors ZEB1 and ZEB2, Genes & Development 22(7) (2008), 894–907.

37.

Vandenboom

T.G.

Kong

et al., Up-regulation of miR-200 and let-7 by natural agents leads to the reversal of epithelial-to-mesenchymal transition in gemcitabine-resistant pancreatic cancer cells, Cancer Research 69(16) (2009), 6704–6712.

38.

Ohuchida

Mizumoto

et al., MicroRNA, hsa-miR-200c, is an independent prognostic factor in pancreatic cancer and its upregulation inhibits pancreatic cancer invasion but increases cell proliferation, Molecular Cancer 9(1) (2010), 169.

39.

Zhu

et al., Serum expression levels of microRNA-382-3p, -598-3p, -1246 and -184 in breast cancer patients, Oncology Letters 12(1) (2016), 269.

40.

Han

Liu

et al., High expression of UBD correlates with epirubicin resistance and indicates poor prognosis in triple-negative breast cancer, Oncotargets & Therapy 8 (2015), 1643–1649.

41.

Eckelpassow

J.E.

Serie

D.J.

Bot

B.M.

et al., ANKS1B is a smoking-related molecular alteration in clear cell renal cell carcinoma, BMC Urology 14(1) (2014), 14.

42.

Quirk

J.T.

Natarajan

et al., Cigarette smoking and the risk of bladder cancer in men and women, Tobacco Induced Diseases 2(3) (2004), 141.

43.

Schmidtke

Kalveram

Weber

et al., The UBA domains of NUB1L are required for binding but not for accelerated degradation of the ubiquitin-like modifier FAT10, Journal of Biological Chemistry 281(29) (2006), 20045–20054.

44.

Abdelzaher

and Kotb

A.F.

, High coexpression of runt-related transcription factor 2 (RUNX2) and p53 independently predicts early tumor recurrence in bladder urothelial carcinoma patients, Applied Immunohistochemistry & Molecular Morphology Aimm 24(5) (2015), 345.

45.

Ehrig

Leivo

Argraves

W.S.

et al., Merosin, a tissue-specific basement membrane protein, is a laminin-like protein, Proceedings of the National Academy of Sciences of the United States of America 87(9) (1990), 3264–3268.

46.

Ying

S.Y.

, Inhibins and activins: Chemical properties and biological activity, Experimental Biology and Medicine 186(3) (1987), 253–264.

Clinical and RNA expression integrated signature for urothelial bladder cancer prognosis

Abstract

BACKGROUND:

PURPOSE:

MATERIALS AND METHODS:

RESULTS:

CONCLUSIONS:

Keywords

1. Introduction

2. Materials and methods

2.1 Patient characteristics and integrated RNA profiles

2.2 Prognostic clinical factors

2.3 Prognostic mRNAs and miRNAs

2.4 Prognosis index and integrated clinical-RNA signature

2.5 Model evaluation

2.6 GO and pathway of mRNA and targets of miRNAs

Table 1 Univariate survival analysis and multivariate Cox regression of demographic and clinical variables

3.1 Clinical characteristics correlated with UBC patient survival

Table 3 Multivariate Cox regression of demographic and clinical variables

3.4 Network analysis

4. Discussion

Footnotes

Acknowledgments

Conflict of interest

Supplement material

References

Table 1
Univariate survival analysis and multivariate Cox regression of demographic and clinical variables

Table 3
Multivariate Cox regression of demographic and clinical variables