Expressing mRNA profiles and related clinical data of human BCa were downloaded from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/) (Barrett et al., 2013). Dataset GSE32894 performed on Illumina Human HT-12 V3.0 expression bead chip was used as the training set (Sjödahl et al., 2012). Dataset GSE13507 performed on Illumina human-6 v2.0 expression bead chip (Kim et al., 2010; Lee et al., 2010) and mRNA expression profiles of BLCA patients were obtained from the TCGA data portal (https://gdc-portal.nci.nih.gov/) (Ye et al., 2019) and were used as another test set. Prognostic data for all TCGA survival analyses were obtained from published papers (Liu et al., 2018).

We used RMA background correction for the raw expression data for the microarray analyses at first, and log₂ transformation and normalization were employed for processed signals. Then, we used the “affy” R package to summarize the median-polish probe sets. The Affymetrix annotation files annotated probes. For TCGA BLCA data, the gene expression data were based on the RNA-sequencing technology of IlluminaHiseq.

Firstly, we excluded samples without exact survival data. By applying the univariate hazard Cox regression analysis with survival as a dependent characteristic, the correlation between each gene expression profile and cancer-specific survival in patients was evaluated based on the training dataset (GSE32894). Here, we identified genes with p < 1E-6 of cancer-specific survival as prognostic gene signatures and then performed LASSO regression analysis. Genes selected from LASSO regression analysis were taken as the candidate factors, and then were subjected to perform multivariate hazard Cox regression analysis in the training dataset with cancer-specific survival as the dependent prognostic influence factor. The risk score was developed based on a linear combination of the mRNA expression level weighted by the estimated regression coefficient generated from the multivariate hazard Cox regression analysis. The formula of risk score for each patient was calculated as follows: Risk score = βgene1 × exprgene1 + βgene2 × exprgene2+ ··· + βgeneN × exprgeneN, in which N is the number of prognostic gene signatures, expr represents the expression profiles of gene signatures, and β means the estimated regression coefficient of gene signatures derived from the multivariate hazard Cox regression analysis. Then, the gene signatures could calculate a risk score for each patient, and we could divide the patients into two (high- and low-risk) groups according to the median risk score. The Kaplan–Meier analysis was used to evaluate the cancer-specific survival distributions by the R “survival” package. Then, another two independent datasets were used to perform the test of the prognostic signature. GSE13507 was used to test the cancer-specific survival and TCGA BLCA data were used to test the disease-specific survival distribution. Moreover, we performed univariate Cox regression and multivariate Cox regression analysis to further verify the prognostic model’s accuracy and precision by integrating clinical features (including gender, age, tumor stage, tumor grade, and progression).

The “limma” R package was utilized to screen the distinguishingly expressed genes between high-risk and low-risk patients. The SAM (significance analysis of microarrays) with FDR (false discovery rate) < 0.05 and |log₂ fold change (FC)| > 1 were set as the cutoff, and the DEGs were applied to further analysis.

Gene ontology (GO) analysis (here, we chose the biological process) was accomplished using the R package cluster Profiler to observe the potential functions of DEGs. p < 0.05 was set as the cutoff criterion.

To further analyze the potential function, the training set was performed into two groups according to the median risk score. For use with GSEA software (https://www.gsea-msigdb.org/gsea/index.jsp) (Subramanian et al., 2005), the collection of annotated gene sets of h.all.v6.1.symbol.gmat [Hallmarks] in Molecular Signatures Database (MSigDB, http://software.broadinstitute.org/gsea/msigdb/index.jsp) was chosen as the reference gene sets (Subramanian et al., 2005; Croken et al., 2014). We selected the gene sets enriched in high-risk groups or high expression level groups, and p < 0.05 was chosen as the cutoff criteria.

To further evaluate the gene expression level between normal bladder and BCa tissues, we used an online database GEPIA2 (http://gepia2.cancer-pku.cn/) (Tang et al., 2019). Moreover, the test set GSE13507 was used to compare the differences between normal bladder mucosae, bladder mucosae surrounding cancer, primary non-muscle invasive BCa, primary muscle invasive BCa, and recurrent non-muscle invasive tumor.

Total RNA was extracted from the nontumorous immortalized bladder cell line (SV-HUV-1) and four BCa cell lines (5637, T24, UM-UC3, and J82) using HiPure Total RNA Mini Kit (Cat. #R4111-03, Magen, China) according to the manufacturer’s instruction. The reverse transcription process was carried out with the ReverTra Ace qPCR RT Kit (Toyobo, China). The expressions of six genes were normalized to GAPDH expression. The primer sequences are listed as Supplementary Table S1.

Univariate hazard Cox regression, LASSO regression, and multivariate hazard Cox regression analyses were performed to identify the prognostic factors and to establish a prognostic model. The survival curve was drawn by the Kaplan–Meier method and compared by log-rank test. ROC curve was used to evaluate the predictive power of the prognostic index. Univariate Cox regression analysis and multivariate Cox regression analysis were performed to further verify the independent prognostic value of the prognostic signature. The statistical significance of differences in qRT-PCR was compared using the Student’s t-test as appropriate. Bioinformatic analysis was done in the R language (version 3.6.2) and p < 0.05 was considered as statistically significant at two sides.

The flow chart of recognition and validation of the six-gene signature is shown as Figure 1. Originally, we employed the univariate hazard Cox regression analysis to assess the connection between all gene expressions and patients’ cancer-specific survival in the training dataset (GSE32894) (Figure 1). Moreover, the result revealed that there were 60 genes significantly associated with prognosis (p < 1E-6), which were defined as prognostic genes. Then, the candidate genes were performed by LASSO regression (Figures 2A,B), and CDK4, GUCY1A2, NMMT, E2F7, ZNF415, HTR2A, NUAK1, COL11A1, THOP1, TNFRSF6B, BCAT1, CBX2, CTRC, DHRS2, BDKRB2, YIF1B, and SLC22A16 were screened. Among these prognostic genes, only three genes (ZNF415, HTR2A, and DHRS2) with higher expression were correlated with more prolonged survival [whose z (coefficient) < 0], whereas other genes (CDK4, GUCY1A2, NMMT, E2F7, NUAK1, COL11A1, THOP1, TNFRSF6B, BCAT1, CBX2, CTRC, BDKRB2, YIF1B, and SLC22A16) with higher expression were lined with shorter survival [whose z (coefficient) > 0].

Multivariate hazard Cox regression analysis was further used to analyze those 17 prognostic genes and then selected genes independently related to cancer-specific survival. Eventually, we screened six genes (CDK4, E2F7, COL11A1, BDKRB2, YIF1B, and ZNF415) as the independent factor and established a prognostic model (Figure 2C). Via integrating the expression of those six genes and the estimated regression coefficient, we then obtained the following calculation model: Risk score = (1.43215589574675 × expression of CDK4) + (0.921330280956022 × expression of E2F7) + (−1.04548254381182 × expression of ZNF415) + (0.814780461026126 × expression of COL11A1) + (0.973314699914422 × expression value of BDKRB2) + (0.964685342668668 × expression value of YIF1B). With the six-gene signature, the risk score for each patient with BCa in the training dataset could be calculated and ranked from the largest to the smallest. Based on the median risk score (0.630561), 224 BCa patients in the training dataset were divided into a high-risk group (n = 112) and a low-risk group (n = 112). There was an obvious difference (p = 6.6956E-08) in patients’ cancer-specific survival between the high-risk and the low-risk groups (Figure 3A). Moreover, we could observe that those ranked into the high-risk group had remarkably shorter survival (median 28.84 months) than those in the low-risk group (median 44.28 months). The time-dependent ROC curve was carried out for 3- and 5-year cancer-specific survival to evaluate the efficacy of the six-gene signature for predicting the cancer-specific survival. The AUCs for the six-gene signature at the cancer-specific survival of 3 and 5 years were 0.96 and 0.967, respectively (Figure 3D). The distribution of the risk score, cancer-specific survival time, and six genes’ expression profiles in the training dataset are shown in Figure 3G, ranked with the increasing risk score. We could find that high-risk patients lived shorter than low-risk patients, and meanwhile, the expression level of patients had a similar trend in five genes (CDK4, E2F7, COL11A1, BDKRB2, and YIF1B), elevating with the increasing risk score, while ZNF415 demonstrated the opposite trend.

Cancer-specific survival of GSE13507 was utilized to test and verify the prognostic efficacy of the six-gene signature for cancer-specific survival prediction; 165 patients of the test set (GSE13507) were classified into the high-risk group (n = 83) and low-risk group (n = 82) according to the same formula generating from GSE32894. The result showed a significant difference (p = 0.0080; median 29.37 vs. 46.835 months) in cancer-specific survival between high-risk and low-risk groups (Figure 3B). The AUC for the six-gene signature was 0.744 and 0.748 at the cancer-specific survival of 3 and 5 years, respectively, in the test set (GSE13507) (Figure 3E). The distribution of the risk score, cancer-specific survival time, and six genes’ expression profiles in the test set of GSE13507 are shown in Figure 3H, ranked with the increasing risk score. In addition, the disease-specific survival of TCGA was used to verify the accuracy of the six-gene signature. As shown in Figure 3C, patients in high risk had a lower survival rate than those in low risk (p = 0.0041). The AUC for the six-gene signature was 0.576 and 0.606 at the disease-specific survival of 3 and 5 years, respectively (Figure 3F). The distribution of the risk score, disease-specific survival time, and six genes’ expression profiles in TCGA are shown in Figure 3I. Above all, the results indicated the good reliability and reproducibility of the six-gene prognostic model for forecasting cancer-specific survival for patients with BCa.

In order to explore whether the prognostic index is an independent prognostic factor, we conducted univariate Cox regression analysis and multivariate Cox regression analysis by integrating several clinicopathological characteristics, including gender, age, tumor stage, tumor grade, and progression. The results indicated that prognostic signature was significantly associated with the cancer-specific survival of BLCA not only in univariate Cox regression analysis (p < 0.001) (Figure 4A), but also in multivariate Cox regression analysis (p < 0.001) (Figure 4B). In summary, the six-gene prognostic model can be seen as an independent prognostic indicator of BLCA.

Subsequently, the correlation of the six-gene signature with clinicopathological features and its prognostic significance were analyzed in the training set and two test sets. We observed that the signature was significantly correlated with BCa divided by T-stage in GSE32894 and GSE13507 (Figures 5A,D) grade in all sets (Figures 5B,E,H). In addition, we found that it was also associated with molecular subtype in GSE32894 (Figure 5C), pathological stage in TCGA (Figure 5G) and progression in test sets GSE13507 (Figure 5F) and TCGA (Figure 5I).

Furthermore, to assess the prognostic value of the six-gene index, the stratified analyses were performed by using clinical information including age, gender, tumor grade, tumor stage, node status, and tumor progression. All 224 BCa patients were firstly stratified by age into the younger dataset (<65 years old, n = 70) and the elder dataset (≥65 years old, n = 154), by gender into a female dataset (n = 61) and male dataset (n = 163), and by tumor grade into grade 1–2 (n = 129) and grade 3 (n = 93). The prognostic power of the six-gene signature was significant in the younger dataset, the elder dataset, the female dataset, the male dataset, the grade 1–2 dataset, and the grade 3 dataset (Figures 6A–F). Based on the tumor stage, patients were categorized into low stage (Ta and T1, n = 173) and high stage (T2–T4, n = 51). Meanwhile, patients were also stratified by node status into N0 (n = 26) and N+ (n = 20) and by tumor progression status into non-tumor progression dataset (n = 211) and tumor progression dataset (n = 13). Interestingly, a similar significant prognostic value could be observed in the high-stage dataset and patients without progression dataset (Figures 6H,I). Otherwise, the prognosis of the low-stage dataset, and N0 and N+ and tumor progression datasets had no significance (data not shown), which may be due to the limited patients.

To investigate the potential function of the six prognostic genes, samples in the training set GSE32894 were divided into two groups according to the risk score. Under the threshold of FDR < 0.05 and |log₂ FC| > 1, a total of 82 DEGs were screened (54 downregulated and 28 upregulated). The volcano plot presented the differential expressed signatures between high- and low-risk groups (Figure 7A).

The biological process of down- and upregulated genes between high- and low-risk groups is visualized in Figures 6B,C, respectively. In the low-risk group, the biological process was enriched in lipid catabolic process, lipid transport, lipid localization, steroid metabolic process, regulation of macrophage-derived foam cell differentiation, triglyceride catabolic process, macrophage-derived foam cell differentiation, foam cell differentiation, neutral lipid catabolic process, and acylglycerol catabolic process. In the high-risk group, the biological processes were significantly enriched in the nuclear division, organelle fission, mitotic nuclear division, chromosome segregation, sister chromatid segregation, nuclear chromosome segregation, mitotic sister chromatid segregation, regulation of mitotic nuclear division, regulation of nuclear division, and regulation of chromosome segregation. Moreover, GSEA analysis was performed, and it revealed that high-risk samples were enriched in G2/M checkpoint, E2F targets, mitotic spindle, mTOR signaling, spermatogenesis, EMT, DNA repair, PI3K/AKT/mTOR signaling, UPR, and MYC targets V2 (Figure 7D).

The results of qRT-PCR and expression profiles between the normal bladder and BCa tissues of six signatures are shown in Figure 7. Compared with normal bladder epithelial cell line (SV-HUV-1), the level of CDK4, E2F7, COL11A1, BDKRB2, and YIF1B was upregulated in most BCa cell lines (Figures 8A,C,E,G,I). On the contrary, the level of ZNF415 (Figure 7K) was downregulated, compared with SV-HUV-1, which were in line with our above contents. In GEPIA2, the expression of CDK4, E2F7, COL11A1, and YIF1B was upregulated in BCa tissues compared with normal bladder tissues (Figures 8B,D,F,J), while the level of BDKRB2 and ZNF415 showed an opposite outcome (Figures 7H,L). FDR < 0.05 and |log₂ FC| > 1 were used as thresholds for judging the significance of gene expression differences in GEPIA2. The results of qRT-PCR were roughly in line with the consequences of GEPIA2 and our above contents that higher expression was related with shorter survival, such as CDK4 and E2F7, and that higher expression was connected with longer survival, for instance, ZNF415.