Identification and validation of an individualized autophagy-clinical prognostic index in bladder cancer patients.

Purpose: Autophagy is a major catabolic system by which eukaryotic cells undergo self-degradation of damaged, defective, or unwanted intracellular components. An abnormal autophagic level is implicated in the pathogenesis of multiple diseases, including cancers. The aim of this study is to explore the prognostic value of autophagy in bladder cancer (BC), which is a major cause of cancer-related death globally. Patients and methods: First, 27 differentially expressed autophagy-related genes (ARGs) were identified in BC patients based on The Cancer Genome Atlas (TCGA) database. Functional enrichment analyses hinted that autophagy may act in a tumor-suppressive role in the initiation of BC. Then, the Cox proportional hazard regression model were employed to identify three key prognostic ARGs (JUN, MYC, and ITGA3), which were related with overall survival (OS) significantly in BC. The three genes represented important clinical significance and prognostic value in BC. Then a prognostic index (PI) was constructed.
Results: The PI was constructed based on the three genes, and significantly stratified BC patients into high- and low-risk groups in terms of OS (HR=1.610, 95% CI=1.200-2.160, P=0.002). PI remained as an independent prognostic factor in multivariate analyses (HR=2.355, 95% CI=1.483-3.739, P<0.001). When integrated with clinical characteristics of age and stage, an autophagy-clinical prognostic index (ACPI) was finally validated, which had improved performance in predicting OS of BC patients (HR=2.669, 95% CI=1.986-3.587, P<0.001). The ACPI was confirmed in datasets of GSE13507 (HR=7.389, 95% CI=3.645-14.980, P<0.001) and GSE31684 (HR=1.665, 95% CI=0.872-3.179, P=0.122).
Conclusion: This study provides a potential prognostic signature for predicting prognosis of BC patients and molecular insights of autophagy in BC.

Introduction

Autophagy, also known as type II programmed cell death, is a major catabolic system by which eukaryotic cells undergo self-degradation of damaged, defective, or unwanted intracellular components.1,2 This is, in part, to quality control of intracellular organelles by continually renewing fresh, better-quality ones. Therefore, stability of cellular renovation, homeostasis, and maintaining physiological level are inseparable from autophagy. An abnormal autophagic level implicated in the pathogenesis of multiple diseases, including inflammation, neurodegenerative diseases, and tumors.3–7 However, the knowledge of autophagy-related mechanism in cancer is still rudimentary and inconclusive. Due to the complex function of autophagy in cancer, the further research on the relation of autophagy and tumors, underlying biological process, and then to apply this knowledge in well-designed therapeutic strategy could be valuable in the new route of cancer therapy. Even whether autophagy is a friend or a foe for cancers cannot draw reliable conclusions for now.8–10

Bladder cancer (BC) is a major cause of cancer-related death globally, causing 165,100 deaths per year.11 In the United States, thee were an estimated 79,030 newly-diagnosed cases in 2017 and 16,870 patients who succumbed to BC.12 Recently, several studies reported that autophagy could be an indispensable mechanism of the onset and progression of BC, which provided a new direction for the clinical management of BC.13–17 Su et al18 observed increased autophagic proteins in high grade urothelial bladder carcinoma, which were regulated via AMPK activation and mTOR inhibition for tumor cells survival, and inhibition of autophagy led to cancer cell death. Some studies have also suggested that targeting autophagy could improve sensitivity to anti-bladder cancer chemotherapy agents.14,15 Thus, exploring the appropriate molecular biomarkers focused on autophagy has attractive value in estimating the deterioration of BC reliably, and may be an important means of fighting BC.

Here we examined the correlation between expression profiles of autophagy-related genes (ARGs) and clinical outcome in 412 BC patients from The Cancer Genome Atlas (TCGA) and developed prognostic index (PI) as an independent index for overall survival (OS) prognosis based on ARGs. To leverage the complementary value of molecular and clinical characteristics, we integrated the PI with clinical factors to build a composite autophagy-clinical prognostic index (ACPI), which allowed us to improve the prognostic efficiency of BC patients. Further validation based on other databases evidently support our risk score model. These findings could also provide an effective multi-dimensional biomarker strategy that would be effective in monitoring autophagy and predicting the prognosis in BC patients.

Materials and methods

Data acquisition

The Human Autophagy Database (HADb, http://www.autophagy.lu/index.html) is an autophagy-dedicated database aiming to reserve human genes involved in autophagy. A variety of ARGs were obtained from the database. RNA-sequencing (RNA-seq) data of ARGs and the clinical information of the bladder urothelial cancer (BLCA) cohort were downloaded and extracted from the TCGA data portal.

Differentially expressed ARGs enrichment analysis

EdgeR package in R statistical software was applied to estimate differentially expressed ARGs between BC and non-tumor samples. Genes exhibiting at least 2-fold changes corresponding to an adjusted P-value less than 0.05 were selected as the significantly differentially expressed ARGs. Then, we performed a series of gene functional enrichment analyses to find the major biological attributes of these genes, including gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG). The Database for Annotation, Visualization, and Integrated Discovery (DAVID, https://david.ncifcrf.gov/), a widely used functional annotation tool, was used to identify enriched GO and KEGG themes. To provide high-dimensional information, the GOplot package of R was performed to concentrate on the visualization of enrichment terms.

Construction of an individualized prognostic index based on ARGs

ARGs expression profiles downloaded from TCGA were normalized by [log2(count+1)] transformed. Univariate Cox regression analyses were performed to select the ARGs whose expression profiles were significantly associated with BC patients’ OS. Subsequently, these survival-related genes were subjected to a multivariate Cox regression analysis to remove the genes that might not be an independent indicator in prognosis monitoring. Finally, several prognostic ARGs were obtained and the PI composed of these genes was developed. The formula of PI based on a linear combination of the relative expression level of genes multiplied regression coefficients, which represented the relative weight of genes in the multiple Cox analysis. BC patients were separated into high- and low-risk groups by the median PI value as the risk cutoff value. The survival curves were plotted by Kaplan–Meier (K–M) method, and differences in the survival rates between high- and low-risk groups were assessed using the log-rank test.

To investigate if the autophagy-related PI could be an independent predictor of OS in the TCGA cohort of BC patients, the multivariate Cox regression analysis was conducted. The PI, age, gender, tumor subtype, pathological stages, and histological grades were used as covariates. Age, stage, and PI were coded as continuous variables. Specifically, stage was coded as I=1, II=2, III=3, and IV=4. The risk factors of gender, subtype, and histologic grade are male, non-Papillary, and high grade.

Statistical analysis

All statistical analyses were conducted using SPSS 24.0 (Chicago, IL, USA) and R 3.3.1 (https://www.r-project.org/). R, GraphPad Prism 5 (San Diego, CA, USA), and OriginPro 2017 (Northampton, MA, USA) were performed to draw plots. Univariate Cox regression analyses were used to evaluate the association between expression profiles and OS. The Multivariate Cox proportional hazards regression model was used to construct the PI and ACPI model based on the factor correlated with survival. Receiver operating characteristic (ROC) curve and the corresponding area under the ROC curve (AUC) for each dataset to measure the prognostic value of ACPI were performed by the package of “survivalROC” in R. All statistical significance was defined as a P-value less than 0.05. Meta-analysis of the selected genes for ACPI was performed based on the expression data of BC and non-BC tissues from Oncomine (https://www.oncomine.org/) with softwares of STATA (version 12.0) and Meta-DiSc (version 1.4).

Results

Differentially expressed ARGs

Altogether RNA-seq and clinical data of 414 BC tissue samples and 19 non-tumor samples were downloaded from TCGA. Among these patients, a total of 408 primary BC patients with gene expression data and clinical follow-up information was involved in the current study. Expression values of 234 ARGs were extracted. Considered as the criteria of a FDR <0.05 and |log2(Fold Change)|>1, we finally obtained nine up-regulated and 18 down-regulated ARGs (Figure 1). Furthermore, scatter plots were visualized to display the expression pattern of the 27 differentially expressed ARGs between BC and non-tumor tissues (Figure 2). Scatter plots displayed expression patterns of 18 down-regulated genes (FOS, JUN, HSPB8, CDKN1A, ITPR1, TP53INP2, PPP1R15A, DLC1, BAG3, MYC, GABARAPL1, BLC2, CCL2, PRKN, NAMPT, CXCR4, NRG2, and CX3CL1) and nine up-regulated genes (BIRC5, CDKN2A, BID, EVA1A, TP73, RGS19, EIF4EBP1, ITGB4, and ITGA3).

Figure 1

Differentially expressed autophagy-related genes (ARGs) between bladder cancer (BC) and normal bladder tissues. (A) The volcano plot for the 234 ARGs from the TCGA data portal. Red indicates high expression and green low expression. Blue shows those genes showed no difference between BC and normal bladder tissues. (B) Hierarchical clustering of differentially expressed ARGs expression levels.

Abbreviations: FDR, false discovery rate; TCGA, The Cancer Genome Atlas.

Figure 2

The expression patterns of 27 autophagy-related genes (ARGs) in bladder cancer types and paired non-tumor samples. Each red dot represents a distinct tumor sample and blue a non-tumor sample. The red bar above the gene name shows a significantly high expression and the blue bar a low expression.

Functional annotation of the differentially expressed ARGs

Functional enrichment analysis of the 27 differentially expressed ARGs offered that the biological understanding of these genes. The GO terms function and KEGG pathway enrichment of these genes were summarized in Table 1. According to the results of DAVID, we found that the top enriched GO terms for biological processes were: response to drug, response to gamma radiation, and apoptotic process; and for cellular components were: cytosol, protein complex, and mitochondrion. On the basis of molecular function, genes were mostly enriched in terms of transcription factor binding, ubiquitin protein ligase binding, and protein heterodimerization activity. The overview schematic of the analysis results is displayed in Figure 3. Besides, in the KEGG pathway enrichment analysis for the differentially expressed ARGs, these genes were shown to be notably associated with Pathways in cancer, colorectal cancer, Hepatitis B, ErbB signaling pathway, p53 signaling pathway, and so on. As shown in Figure 4A, the Z-score of enriched pathways less than zero indicated that most of the cancer pathways were more likely to be decreased. The heatmap of the relationship between ARGs and pathways was also displayed (Figure 4B).

Table 1

GO and KEGG analysis of differentially expressed autophagy-related genes

Category	ID	Term	P-value	Genes
Biological Process	GO:0042493	Response to drug	4.44E−06	FOS, CDKN1A, JUN, BCL2, ITGA3, MYC, TP73
Biological Process	GO:0010332	Response to gamma radiation	1.27E−05	CCL2, BCL2, MYC, TP73
Biological Process	GO:0006915	Apoptotic process	1.36E−05	DLC1, EVA1A, CDKN2A, CXCR4, BCL2, BIRC5, PPP1R15A, TP73
Biological Process	GO:0007346	Regulation of mitotic cell cycle	2.78E−05	CDKN1A, BIRC5, MYC, TP73
Biological Process	GO:0043524	Negative regulation of neuron apoptotic process	4.06E−05	CCL2, JUN, BCL2, BIRC5, TP73
Biological Process	GO:0007050	Cell cycle arrest	5.25E−05	CDKN1A, CDKN2A, PPP1R15A, MYC, TP73
Biological Process	GO:0006974	Cellular response to DNA damage stimulus	2.36E−04	CDKN1A, BCL2, PPP1R15A, MYC, TP73
Biological Process	GO:0001836	Release of cytochrome c from mitochondria	5.28E−04	BID, JUN, BCL2
Biological Process	GO:0045893	Positive regulation of transcription, DNA-templated	8.50E−04	FOS, CDKN2A, JUN, MYC, TP53INP2, TP73
Biological Process	GO:1900740	Positive regulation of protein insertion into mitochondrial membrane involved in apoptotic signaling pathway	9.02E−04	BID, BCL2, TP73
Cellular Component	GO:0005829	Cytosol	3.94E−06	BID, DLC1, NAMPT, GABARAPL1, BIRC5, TP73, FOS, CDKN1A, EIF4EBP1, CDKN2A, JUN, BCL2, BAG3, MYC, PPP1R15A, TP53INP2
Cellular Component	GO:0043234	Protein complex	0.002235	EIF4EBP1, CDKN1A, CDKN2A, MYC, ITPR1
Cellular Component	GO:0005739	Mitochondrion	0.007959	BID, GABARAPL1, CDKN2A, BCL2, PPP1R15A, MYC, TP73
Cellular Component	GO:0005654	Nucleoplasm	0.008942	NAMPT, FOS, EIF4EBP1, CDKN1A, CDKN2A, JUN, HSPB8, BIRC5, MYC, TP73
Cellular Component	GO:0005741	Mitochondrial outer membrane	0.017612	BID, BCL2, PPP1R15A
Cellular Component	GO:0030054	Cell junction	0.024203	NAMPT, CXCR4, ITGB4, TP73
Cellular Component	GO:0005783	Endoplasmic reticulum	0.025005	FOS, GABARAPL1, BCL2, PPP1R15A, ITPR1
Cellular Component	GO:0005667	Transcription factor complex	0.028532	FOS, JUN, TP73
Cellular Component	GO:0005794	Golgi apparatus	0.028578	GABARAPL1, HSPB8, RGS19, PPP1R15A
Cellular Component	GO:0008305	Integrin complex	0.036412	ITGB4, ITGA3
Cellular Component	GO:0005515	Protein binding	3.97E−05	BID, DLC1, NAMPT, GABARAPL1, ITGB4, RGS19, BIRC5, ITGA3, CX3CL1, TP73, ITPR1, FOS, CDKN1A, EIF4EBP1, CDKN2A, CXCR4, JUN, HSPB8, BCL2, BAG3, MYC, PPP1R15A, TP53INP2
Molecular Function	GO:0008134	Transcription factor binding	4.25E−05	FOS, CDKN2A, JUN, BCL2, MYC, TP73
Molecular Function	GO:0031625	Ubiquitin protein ligase binding	6.65E−04	BID, GABARAPL1, CDKN1A, CXCR4, BCL2
Molecular Function	GO:0046982	Protein heterodimerization activity	0.003898	FOS, JUN, BCL2, BIRC5, ITGA3
Molecular Function	GO:0001077	Transcriptional activator activity, RNA polymerase II core promoter proximal region sequence-specific binding	0.004395	FOS, JUN, MYC, TP73
Molecular Function	GO:0043565	Sequence-specific DNA binding	0.005716	FOS, JUN, BCL2, MYC, TP73
Molecular Function	GO:0000978	RNA polymerase II core promoter proximal region sequence-specific DNA binding	0.013453	FOS, JUN, MYC, TP73
Molecular Function	GO:0004861	Cyclin-dependent protein serine/threonine kinase inhibitor activity	0.015533	CDKN1A, CDKN2A
Molecular Function	GO:0019901	Protein kinase binding	0.015684	CDKN2A, RGS19, PPP1R15A, TP73
Molecular Function	GO:0042802	Identical protein binding	0.020127	JUN, HSPB8, BCL2, BIRC5, TP73
KEGG PATHWAY	hsa05200	Pathways in cancer	3.16E−07	BID, FOS, CDKN1A, CDKN2A, CXCR4, JUN, BCL2, BIRC5, ITGA3, MYC
KEGG PATHWAY	hsa05210	Colorectal cancer	2.06E−05	FOS, JUN, BCL2, BIRC5, MYC
KEGG PATHWAY	hsa05161	Hepatitis B	3.48E−05	FOS, CDKN1A, JUN, BCL2, BIRC5, MYC
KEGG PATHWAY	hsa04012	ErbB signaling pathway	7.87E−05	EIF4EBP1, CDKN1A, JUN, NRG2, MYC
KEGG PATHWAY	hsa04115	p53 signaling pathway	7.56E−04	BID, CDKN1A, CDKN2A, TP73
KEGG PATHWAY	hsa04151	PI3K-Akt signaling pathway	0.001961	EIF4EBP1, CDKN1A, BCL2, ITGB4, ITGA3, MYC
KEGG PATHWAY	hsa04668	TNF signaling pathway	0.002843	FOS, CCL2, JUN, CX3CL1
KEGG PATHWAY	hsa05166	HTLV-I infection	0.004594	FOS, CDKN1A, CDKN2A, JUN, MYC
KEGG PATHWAY	hsa05219	Bladder cancer	0.00551	CDKN1A, CDKN2A, MYC
KEGG PATHWAY	hsa04921	Oxytocin signaling pathway	0.008685	FOS, CDKN1A, JUN, ITPR1

Abbreviations: GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Figure 3

The bubble plot of enriched gene ontology (GO) terms. The z-score is assigned to the x-axis, and the negative logarithm of the P-value to the y-axis, as in the barplot (the higher the more significant). The size of the displayed circles is proportional to the number of genes assigned to the term. Greed circles correspond to the biological process, red indicates the cellular component, and blue shows the molecular function category.

Figure 4

Kyoto Encyclopedia of Genes and Genomes analysis of differentially expressed autophagy-related genes (ARGs). (A) The outer circle shows a scatter plot for each term of the logFC of the assigned genes. Red circles display up-regulation, and blue ones down-regulation. (B) The heatmap of the relationship between ARGs and pathways. The color of each block depends on the logFC values.

Identification of prognostic ARGs

The relationships between the expression profiles of 27 differentially expressed ARGs and OS were assessed based on the data obtained from TCGA, resulting in four prognosis-related ARGs. In order to improve the robustness, four prognosis-related ARGs (JUN, MYC, ITGA3, and NAMPT) were selected for further multivariate Cox regression model by SPSS 24.0 (Table 2). However, the gene of NAMPT showed no significant prognostic value with P>0.05. Finally, three genes including JUN, MYC, and ITGA3 were identified to develop the PI model (Table 2). The results from K-M analysis indicated that the up-regulation of JUN was strongly correlated with the inferior OS of BC patients (HR=1.925, 95% CI=1.325–2.798, P<0.001; Figure 5A). Also, MYC overexpression lead to worse OS (HR=1.931, 95% CI=1.426–2.614, P<0.001; Figure 5C). On the contrary, up-regulated ITGA3 indicated BC patients has a longer survival time (HR=0.659, 95% CI=0.487–0.893, P=0.007; Figure 5E). According to the median value of the three genes to group, the same trend was obtained (Figures 5B, D, and F).

Table 2

Expression and Cox regression analysis data of the prognosis-related ARGs in bladder cancer by TCGA

Gene		Expression		Univariate Cox			Multivariate Cox
mean±SD	P-value	HR	z	P-value	HR	z	P-value
JUN	Tumor	13.002773±1.164031	<0.0001	1.198762	2.651598	0.008011	1.1786	2.38	0.017
	Non-tumor	15.027896±1.185059
MYC	Tumor	11.541411±1.554496	<0.0001	1.14059	2.606545	0.009146	1.1682	2.97	0.003
	Non-tumor	13.287664±1.265550
ITGA3	Tumor	13.822128±1.515658	<0.0001	0.907077	−2.0139	0.04402	−0.1505	−2.97	0.003
	Non-tumor	12.975430±1.091134
NAMPT	Tumor	12.330813±1.009040	<0.0001	1.172207	2.001379	0.045352			0.340
	Non-tumor	13.308780±1.459465
FOS	Tumor	13.343147±1.608937	<0.0001	1.063323	1.267733	0.204893
	Non-tumor	16.672991±1.404332
HSPB8	Tumor	10.106659±2.154647	<0.0001	1.068919	1.911278	0.055969
	Non-tumor	13.749171±1.702493
CDKN1A	Tumor	13.124746±1.128535	<0.0001	1.01594	0.22579	0.821364
	Non-tumor	14.389922±1.347070
ITPR1	Tumor	9.5318665±1.169013	<0.0001	0.997873	−0.03415	0.972754
	Non-tumor	11.929274±1.585768
TP53INP2	Tumor	10.462387±0.976958	<0.0001	1.119264	1.518161	0.128974
	Non-tumor	12.653637±1.109542
BIRC5	Tumor	11.387251±1.094559	<0.0001	1.06901	0.897696	0.369347
	Non-tumor	7.9468232±2.565110
PPP1R15A	Tumor	12.441886±0.891543	<0.0001	0.940067	−0.71379	0.475356
	Non-tumor	13.738993±1.001618
DLC1	Tumor	9.5131164±1.257221	<0.0001	1.021996	0.347978	0.727857
	Non-tumor	11.242183±0.995269
BAG3	Tumor	11.719572±0.731981	<0.0001	1.192655	1.702932	0.088581
	Non-tumor	12.704785±1.108736
GABARAPL1	Tumor	11.058632±0.792235	<0.0001	1.051224	0.525272	0.599394
	Non-tumor	12.255963±0.614109
BCL2	Tumor	8.3027604±1.178505	<0.0001	1.051513	0.768008	0.442482
	Non-tumor	10.256868±0.425816
CDKN2A	Tumor	8.9951212±3.453185	<0.0001	0.990449	−0.44257	0.658076
	Non-tumor	6.6235955±1.685343
CCL2	Tumor	9.0438964±1.821760	<0.0001	1.053817	1.273601	0.202805
	Non-tumor	11.704824±2.040782
PRKN	Tumor	5.2978598±1.521700	<0.0001	1.009973	0.201647	0.840193
	Non-tumor	8.1483213±1.213423
BID	Tumor	11.340986±0.737984	<0.0001	0.839371	−1.70367	0.088442
	Non-tumor	10.105586±0.651137
EVA1A	Tumor	6.5437907±2.147278	<0.0001	1.057202	1.528728	0.126332
	Non-tumor	3.8921911±1.267249
CXCR4	Tumor	10.259886±1.380760	<0.0001	1.013263	0.261675	0.793572
	Non-tumor	11.447608±1.958345
NRG2	Tumor	4.7761532±2.145474	<0.0001	0.985697	−0.43265	0.665266
	Non-tumor	7.9643803±1.000720
TP73	Tumor	7.9164301±2.009466	<0.0001	0.984879	−0.3943	0.69336
	Non-tumor	5.9146373±2.363484
RGS19	Tumor	9.9259644±0.729683	<0.0001	1.143028	1.313883	0.188886
	Non-tumor	8.6949825±0.978522
EIF4EBP1	Tumor	11.427944±1.081142	<0.0001	1.046433	0.666103	0.505345
	Non-tumor	9.8783883±0.927159
CX3CL1	Tumor	9.7520038±1.782426	<0.0001	1.034169	0.77976	0.435532
	Non-tumor	12.084799±0.955857
ITGB4	Tumor	14.253622±1.551662	<0.0001	0.936675	−1.33788	0.180934
	Non-tumor	13.188612±1.541228

Abbreviations: ARG, autophagy-related genes; TCGA, The Cancer Genome Atlas.

Figure 5

The correlation between three genes included in prognostic signature and bladder cancer patients’ survival. Kaplan–Meier plots summarize results from analysis of correlation between (A) JUN expression level and patient survival, using best separation, (B) JUN expression level and patient survival, using median separation, (C) MYC expression level and patient survival, using best separation, (D) MYC expression level and patient survival, using median separation, (E) ITGA3 expression level and patient survival, using best separation, (F) ITGA3 expression level and patient survival, using median separation.

The correlations between expression level of the three genes and clinicopathological parameters in BC are summarized in Table 3. We observed significant correlations between JUN overexpression and tumor subtype of non-Papillary (P<0.001), high histological grade (P=0.033), advanced pathologic stage (P<0.001), and advanced pathological T stage (P=0.002). Elevated MYC was closely linked with tumor subtype of non-Papillary (P=0.001) and high histological grade (P=0.005). High expression of ITGA3 occurred in low histological grade (P<0.001), early pathological stage (P=0.033), early pathological T stage (P=0.008), and no lymph node metastasis (P=0.016).

Table 3

Relationships between three key prognostic genes expression and clinical parameters in bladder cancer by TCGA

Parameters		N	JUN expression value			MYC expression value			ITGA3 expression value
M±SD	t value	P-value	M±SD	t value	P-value	M±SD	t value	P-value
Tissues	Tumor	414	13.0028±1.1640	−7.410	<0.001	11.5414±1.5545	−4.822	<0.001	13.8221±1.5157	2.405	0.017
Non-tumor	19	15.0279±1.1851	13.2877±1.2655	12.9754±1.0911
Age	≥60	321	13.0359±1.1327	1.076	0.150	11.6119±1.5304	1.571	0.117	13.7787±1.5272	−1.665	0.097
<60	87	12.8852±1.2530	11.3177±1.6171	14.0830±1.4561
Gender	Male	301	12.9868±1.1674	−0.495	0.621	11.4625±1.5666	−1.898	0.058	13.7850±1.5661	−1.311	0.191
Female	107	13.0515±1.1410	11.7930±1.4901	14.0084±1.3574
Subtype	Non-papillary	271	13.1979±1.0811	4.903	<0.001	11.7343±1.5708	3.405	0.001	13.7762±1.5897	−1.492	0.136
Papillary	132	12.6091±1.2284	11.1787±1.4654	14.0137±1.2921
Tumor status	With tumor	151	13.0248±1.0119	0.459	0.647	11.6697±1.5948	1.351	0.177	13.7276±1.5787	1.763	0.079
Tumor free	201	12.9720±1.2459	11.4571±1.4702	13.9996±1.4314
Histological grade	High grade	384	13.0359±1.1529	2.143	0.033	11.5759±1.5844	3.002	0.005	13.7990±1.5338	−6.996	<0.001
Low grade	21	12.4819±1.1653	11.0637±0.6886	14.8231±0.5669
Pathologic stage	Ⅲ～Ⅳ	274	13.1826±1.0448	4.206	<0.001	11.6576±1.4617	1.944	0.053	13.7353±1.5386	−2.136	0.033
I~II	132	12.6367±1.3030	11.3387±1.7148	14.0770±1.4496
T stage	T3–T4	252	13.1958±1.0249	3.064	0.002	11.6811±1.4423	1.643	0.101	13.7290±1.5227	−2.684	0.008
T1–T2	123	12.8322±1.1819	11.4088±1.6306	14.1726±1.4599
N stage	N1–3	129	13.1086±0.9329	1.265	0.207	11.4309±1.4607	−1.517	0.130	13.6208±1.4944	−2.427	0.016
N0	237	12.9643±1.2176	11.6798±1.5195	14.0205±1.5112
M stage	M1	11	12.8009±1.0016	−0.191	0.849	11.9098±1.8601	1.245	0.215	13.5755±1.3148	−0.792	0.429
M0	196	12.8728±1.2223	11.3199±1.5104	13.9395±1.4919

Abbreviations: M±SD, mean±standard deviation; TCGA, The Cancer Genome Atlas.

Construction and definition of the PI

The formula of PI is as follows: PI=(0.1643 × expression value of JUN)+(0.1555 × expression value of MYC)+(−0.1505 × expression value of ITGA3). It is noticed that the coefficient of ITGA3 is negative, indicating that the expression of JUN and MYC were negatively related with the survival time of BC patients, while the JUN was positively related with OS. Based on the median expression value of PI, the BC patients were stratified into high- and low-risk groups. We also calculated the expression levels of the three prognostic genes between high- and low-risk groups. Remarkably higher expression was noted for JUN and MYC in the high-risk groups, while lower expression was observed for ITGA3 in the high-risk groups (Figure 6). These findings also hint that JUN and MYC were risk factors, while ITGA3 was a protective factor for the progression of BC patients.

Figure 6

Different expression of the three key genes between the high risk group and low risk group.

In the meantime, the relationships between clinicopathological parameters and PI were also investigated. The results of independent sample t-tests showed that the PI values were higher in elder than in younger patients (P=0.009; Figure 7A), higher in non-papillary than in papillary bladder cancer (P<0.001; Figure 7C), higher in TIII–IV than in TI–II (P<0.001; Figure 7D), higher in histological stage III–IV than in I–II (P<0.001; Figure 7G), and higher in high grade than in low grade (P<0.001; Figure 7H). No difference of PI value was observed between male and female (P=0.494; Figure 7B), N1–3 stage and N0 stage (P=0.250; Figure 7E), or M1 and M0 stage (P=0.254; Figure 7F).

Figure 7

The clinicopathological significance of prognostic index (PI) in bladder cancer. PI value in different (A) ages, (B) genders, (C) tumor subtypes, (D) pathological T stages, (E) pathological N stages, (F) pathological M stages, (G) pathological stages, (H) histological grades.

To identify the performance of PI in predicting the clinical outcome of BC patients, the K-M plots were plotted to analyze the different survival time between the high- and low-risk groups. The results of K-M analysis indicated that the median OS for the high-expression group was 734 days; the median OS for the low-expression group was 1,423 days. Patients in the high-risk group suffered significantly worse survival than those in the low-risk group (HR=1.610, 95% CI=1.200–2.160, P=0.002, Figure 8A). Figures 8B–F show the PI distribution of patients in the training dataset, the number of patients in different risk groups, the OS of patients in the TCGA dataset, the number of censor patients, and the heatmap of the three genes expression profiles in the TCGA dataset. Furthermore, PI remained as an independent prognostic indicator for BC patients in multivariate analyses, after adjusting for clinicopathological features such as age, gender, tumor subtype, pathologic stage, and histological grade (HR=2.355, 95% CI=1.483–3.739, P<0.001, Table 4).

Figure 8

Autophagy-related prognostic index (PI) of bladder cancer patients. (A) Kaplan–Meier plot represents that patients in the high-risk group had significantly shorter overall survival time than those in the low-risk group. (B) The PI distribution of patients in the training dataset. (C) The number of patients in different risk groups. (D) The overall survival of patients in the TCGA dataset. (E) The number of censor patients. (F) The heatmap of the three key genes expression profiles in the TCGA dataset.

Abbreviations: TCGA, The Cancer Genome Atlas.

Table 4

Univariate and multivariate analyses of OS in bladder cancer patients of TCGA

Variables	Univariate analysis		Multivariate analysis
Hazard ratio (95%CI)	P-value	Hazard ratio (95% CI)	P-value
Age	1.033 (1.017–1.049)	<0.001	1.029 (1.014–1.045)	<0.001
Gender	0.872 (0.631–1.203)	0.404	0.818 (0.588–1.139)	0.235
Subtype	1.458 (1.030–2.065)	0.033	1.084 (0.756–1.553)	0.661
Pathologic stage	1.707 (1.412–2.065)	<0.001	1.617 (1.321–1.978)	<0.001
Histologic grade	2.968 (0.734–11.995)	0.127	0.931 (0.221–3.918)	0.922
Prognostic index	2.717 (1.764–4.184)	<0.001	2.355 (1.483–3.739)	<0.001

Notes: Age, stage, and prognostic index were coded as continuous variables. Specifically, stage was coded as I=1, II=2, III=3. IV=4. The risk factors of gender, subtype,and histologic grade are male, non-papillary, and high grade.

Abbreviations: OS, overall survival; TCGA, The Cancer Genome Atlas.

Integrated prognostic signature by combining the PI with clinical parameters

Based on the multivariate Cox regression analysis with TCGA dataset, age, stage, and PI were suggested as independent prognostic factors with complementary value. To further improve accuracy of PI in predicting OS of BC patients, we integrated age, pathological stage, and PI to derive an ACPI as (0.028 × age)+(0.467 × stage)+(0.834 × PI score). Similarly, patients were divided into high- and low-risk groups based on the median value of ACPI. As expected, ACPI stratified BC patients into two groups with a significantly different prognosis (HR=2.669, 95% CI=1.986–3.587, P<0.001; Figure 9A). To evaluate how well the ACPI predicts the prognoses of BC patients, the time-dependent ROC curve analysis was carried out. The AUC for the ACPI was 0.689 (Figure 9B), demonstrating the competitive performance of the ACPI for survival prediction in the TCGA dataset. The prognostic value of ACPI was also validated by GSE13507 and GSE31684. Consistent with the findings based on the TCGA dataset, patients in the high-risk group had significantly shorter overall survival than those in the low-risk group based on GSE13507 (HR=7.389, 95% CI=3.645–14.980, P<0.001; Figure 9C), and the AUC for the ACPI was 0.864 (Figure 9D). A similar trend was observed in GSE31684 (HR=1.665, 95% CI=0.872–3.179, P=0.122; Figure 9E), and the AUC for the ACPI was 0.624 (Figure 9F).

Figure 9

The prognostic value of autophagy-clinical prognostic index (ACPI) of bladder cancer patients. (A) Kaplan–Meier (K–M) survival curve showing overall survival (OS) outcomes according to relative high-risk and low-risk patients based on The Cancer Genome Atlas (TCGA) database. (B) Time-dependent ROC curve analysis for survival prediction by the ACPI based on the TCGA database. (C) K–M survival curve showing OS outcomes according to relative high-risk and low-risk patients based on the GSE13507 dataset. (D) Time-dependent ROC curve analysis for survival prediction by the ACPI based on the GSE13507 dataset. (E) K–M survival curve showing OS outcomes according to relative high-risk and low-risk patients based on the GSE31684 dataset. (F) Time-dependent ROC curve analysis for survival prediction by the ACPI based on the GSE31684 dataset.

Meta-analysis

A total of 19 eligible studies were involved, including Blaveri Bladder 2, Modlich Bladder, Sanchez Carbayo Bladder 2, TCGA, GSE3167, GSE13507, GSE76211, GSE2109, GSE7476, GSE30522, GSE31189, GSE37815, GSE52519, GSE65635, GSE37817, GSE100926, GSE24152, GSE19915 (GPL3883 and GPL5186), and GSE40355. The results of meta-analysis and the diagnostic tests of meta-analysis were also updated. The expression of JUN and MYC in BC tissues were lower than that in non-BC tissues (I2JUN=91.6%, PJUN<0.001; I2MYC=90.3%, PMYC<0.001) (Figures 10A and C), while the expression of ITGA3 was opposite (I2ITGA3=92.7%, PITGA3<0.001) (Figure 10E), the same expression tendence with that in TCGA. In addition, the diagnostic tests of meta-analysis showed that the AUC of the sROC of JUN, MYC, and ITGA3 were 0.91, 0.87, and 0.74, respectively (Figures 10B, D, and F). Among the 19 studies involved in meta-analysis, only two microarrays (GSE37137 and GSE35824) and TCGA contained gene expression data from non-muscle invasive bladder cancer (NMIBC) and muscle invasive bladder cancer (MIBC) tissues. Then we performed a meta-analysis to evaluate the expression of JUN, MYC, and ITGA3 between NMIBC and MIBC tissues. The expression of JUN and ITGA3 in MIBC tissues were higher than that in NMIBC tissues (I2JUN=72.1%, PJUN=0.028; I2ITGA3=92.8%, PITGA3<0.001), while the expression of MYC was opposite (I2MYC=67.2%, PITGA3=0.047) (data not shown). The heterogeneity of meta-analysis was significant due to the small size of cases involved in the microarrays and TCGA.

Figure 10

Meta-analysis. (A) Forest plot of JUN expression in bladder cancer with six datasets. (B) sROC curve for JUN expression in bladder cancer with six datasets. (C) Forest plot of MYC expression in bladder cancer with six datasets. (D) sROC curve for MYC expression in bladder cancer with six datasets. (E) Forest plot of ITGA3 expression in bladder cancer with six datasets. (F) SROC curve for ITGA3 expression in bladder cancer with six datasets.

Discussion

BC is a major lethal malignancy worldwide. The stalled advance in molecular targeted therapy and no effective molecular biomarkers for BC prognosis monitoring warrants a better understanding of the molecular mechanisms that underlie this condition.19,20 Exploration of autophagy mechanism opens new perspectives for BC.21–25 However, most research only focused on autophagy via studying a signal gene. To capture the genes necessary for BC from the perspective of autophagy, we screened ARGs and identified three key prognostic ARGs, all of which may offer additional potential therapeutic targets. We further leveraged the complementary value of molecular and clinical characteristics and showed that combining both could provide a more accurate estimation of overall survival in BC. This integrated study of multiple databases contributed to our novel understanding of BC biology and delineated potential therapeutic intervention possibilities.

Given great advances in high-throughput sequencing recently, several large-scale databases emerged, such as TCGA and GEO, which have provided effective measures for selecting gene signatures. In the current study, we deeply mined the expression profiles of ARGs from TCGA and aimed to search molecular biomarkers for detecting the prognosis of BC patients. We first screened 27 differentially expressed ARGs between BC and non-tumor tissues. Considering these genes may be depth involved in the initiation of BC, we performed GO and KEGG analysis of these genes. Interestingly, functional analysis revealed that the most significant KEGG pathway (pathways in cancer) of these enriched genes was decreased. Based on the results, we hypothesized that autophagy may act as the tumor suppressor in the process of tumor initiation. Autophagy caused great concern; of particular interest was its multi-faceted character in cancers. Initially, the tumor-suppressive role of autophagy in cancers was proposed for autophagy inhibited by activation of mutations in oncogenes or inactivation of tumor suppressor genes.26 Furthermore, systemic mosaic deletion of autophagy genes in the setting of certain mouse models can result in the initiation of neoplasia.27 Interestingly, autophagy turns to the guardian of malignant tumor cells after tumors are established.28 However, the role change of autophagy is not immutable and varies in different tumors.29,30

The result of univariate survival analysis revealed that four ARGs were associated with OS in the TCGA database. Further multivariate survival analysis helped us determine three key prognostic ARGs (JUN, MYC, ITGA3) to develop the PI, which could be an independent prognostic indicator for BC patients. JUN encodes c-Jun, which is the first discovered oncogenic transcription factor,31 involving diverse cellular processes, such as cell cycle progression,32 anti-apoptotic, and tumorigenesis. Previous studies have suggested that up-regulation of c-Jun proteins was predictive of inferior OS for BC patients.33 However, well-informed insights of the functional mechanism of JUN in BC only has little coverage. The MYC protein is a multifunctional, nuclear phosphoprotein, and shows its evil face in the progress of a variety of tumors, including BC.34–37 Massari et al38 also found that c-Myc could exert excellent ability in stratifying patients with muscle invasive bladder urothelial carcinoma into high-risk and low-risk groups significantly for survival. In addition, several studies found that c-Myc knockdown could inhibit proliferation, migration, and invasion of bladder cancer cells.37 ITGA3 belongs to a family of the integrins, which triggers cell survival, proliferation, or migration events.39 The present study demonstrated that JUN and MYC overexpression were significantly associated with advanced pathological stage and high grade. Additionally, up-regulation of JUN and MYC indicated inferior OS. The opposite pattern was observed in the relationships between ITGA3 and clinical significance. Hence, we speculated that JUN and MYC may function as major driving forces of tumor progression, while ITGA3 exerted its tumor suppressor role.

To date, some prognostic signature of cancers based on expression profiles were proposed by the aid of advances in a large-scale public database. For example, Bao et al40 analyzed the RNA-Seq data of 234 BC patients from TCGA and managed to obtain a four-lncRNA signature, which exerted a prognosis predicting value. Zhong et al41 also proposed a prognostic signature with six genes as a potential survival prediction marker for ER-positive breast cancer patients. However, these studies only focused on molecular biomarkers and overlooked the traditional clinical parameters. We attached much weight on molecular mechanisms and clinical perspective at once. Thus, the prognostic signature is promising to be converted into clinical application. However, a limitation of this study is its retrospective nature. Due to the lack of enough cases, we failed to evaluate the expression of JUN, MYC, and ITGA3 between NMIBC and MIBC tissues. In addition, other potential prognostic variables correlated to OS in BC, such as body mass index (BMI), residual tumor at tur, neutrophil-to-lymphocyte ratio (NLR), and lymphovascular invasion (LVI), should be investigated. Last, the changes of before and after the treatment, such as chemotherapy or Bacillus Calmette-Guérin (BCG) refractory, should also be considered to find the potential markers for predicting the treatment effect and prognosis.

Conclusion

In conclusion, based on the comprehensive analyses with ARGs expression profiles and corresponding clinical features, three prognostic ARGs (JUN, MYC, and ITGA3) were identified. The genes identified in autophagy pathways also provide new possibilities for bladder cancer therapeutic intervention. By combining molecular signature and clinical characteristics, we constructed a novel risk score model ACPI which can robustly estimate BC patients’ survival. Also, the ACPI risk score model was validated by large sample size. However, further prospective experiments can be expected to test the clinical utility and aid in the search for optimal personalized targeted therapies.