An RNA-sequencing-based transcriptome for a significantly prognostic novel driver signature identification in bladder urothelial carcinoma

Data source and preprocessing

Gene expression profiles (Illumina HiSeq RNA Seq), miRNA mature strand expression profiles (Illumina HiSeq microRNA Seq) level 3 data and phenotype data of bladder urothelial carcinoma (BLCA) (Robertson et al., 2017) from the TCGA database were downloaded through the UCSC Xena portal (https://xena.ucsc.edu/). The mRNA expression profiles and lncRNA expression profiles were derived from the entire gene expression profiles. Genes with missing expression values in > 50% samples were removed. Patients chosen for model constructing met the following criteria: (1) histologic diagnosis of primary BC; and (2) available RNA expression profiles and complete clinic-pathological and follow-up data. After sample filtering, four hundred patients were enrolled for further analysis and divided into 200 samples in the training set and 200 samples in the test set randomly based on the “sample” algorithm. The “sample” algorithm was a basic function in R. The application was generating random samples. This means that sample takes a sample of the specified size from the elements of x using either with or without replacement. We utilized this function to guarantee the randomness of the samples and repeatability of results and avoided the sample-specific results. The selection process of the prognostic signature is shown in Fig. 1.

Preliminary screening of generally changed RNAs

First, RNAs with significantly high variability expression among different patients of primary bladder cancer were preliminary screened. Because these RNAs could be considered as important roles in tumorigenesis. A generally changed gene was selected based on the criteria outlined below:

1. The median expression of gene in every sample was more than 20% of the total median expressions of all RNAs in every sample;

2. The variance of expression of gene in every sample was higher than 20% of the total expressions variances of all RNAs in every sample.

Selection of seed RNAs

All of the selected RNAs with significant changes across patient samples were analysed by univariate Cox regression survival analysis to evaluate the relationship between the expression of each PCG, lncRNA or microRNA and patients’ OS by “survival” package in R in the training set (Grambsch & Therneau, 2000; Therneau, 2015). RNAs found to have significant p values < 0.05 were selected as seed RNAs, indicating prognostic RNAs.

Identification of key RNAs related to prognosis

Considering that a few PCGs, lncRNAs and microRNAs in one model would make the prediction more precise and practical, we selected the top 20 PCGs, lncRNAs, and microRNAs to perform the robust likelihood-based survival models to identify key biomarkers influencing the clinical outcome of cancer by utilizing the “rbsurv” package in R language (Hyung, Cho & Yu, 2009; Kendall, Pollock & Brownie, 1995; Renaud et al., 2015; Wang & Tai, 2009).

Furthermore, we performed the random forests for survival regression and classification (RF-SRC) algorithm by using the “randomForestSRC” package in R language to filter genes until three PCGs, three lncRNAs and three microRNAs were retained to construct a predictive model (Ishwaran & Kogalur, 2010; Ishwaran et al., 2014).

Screening of BLCA subtypes-associated metagenes

Through data preprocessing, we screened out 8969 PCGs, 89 lncRNAs and 582 microRNAs. Before performing NMF, we applied a filtering procedure to retain genes with high variability across 200 patients from the training cohorts because the higher variable genes are informative in the clustering process (Zhao, Zhao & Yan, 2018). The median absolute deviation (MAD) value of each gene was calculated, and we selected the top 500 most variable genes for the clustering analysis. In particular, the 89 lncRNAs expression in the 200 patients were all reserved for the following analysis.

The “NMF” R package (Gaujoux & Seoighe, 2010) was utilized to perform the clustering analysis with the Brunet algorithm. Informative genes were extracted by NMF method. We used NMF to reduce the dimensionality of expression data from thousands of genes to a handful of metagenes (Brunet et al., 2004). The number of clusters k was varied from 2 to 10, and we repeated the clustering process 30 times. The value of k that led to the stable cophenetic correlation coefficient mentioned by Brunet et al. was chosen as the optimal number of clusters. Next, we repeated the clustering 200 times with an optimal k to obtain the associated metagenes and consensus matrix. Subtype-specific genes were singled out using the “extractFeatures” function with the largest row feature scores in the “NMF” package (Gaujoux & Seoighe, 2010).

Identification of BLCA subtypes

We clustered the 200 tumour samples in the training set using the SNF algorithm based on the “CancerSubtypes” R package (Xu et al., 2017). Then, we selected 20 as the number of iterations for the diffusion process. The statistical significance of the clustering (SigClust) was calculated to validate the significance of the clustering results. This calculation assesses the significance of clustering by simulation from a single null Gaussian distribution. The parameter of the “icovest” was set as 3, indicating original background noise threshold estimate. Survival analysis was also conducted to validate the significance and verify the survival patterns between the identified molecular subtypes or in every cancer subtype. We also made a comparison between our classification and the existing subtype classification from the TCGA subtypes (Robertson et al., 2017). The heatmap corresponding to the dendrogram was generated using the heatmap function with SNF algorithm classification, mRNA cluster, lncRNA cluster, microRNA cluster, RPPA cluster, hypermethylation cluster, hypomethylation cluster, mutation process (MSig) cluster, SMG-SCNA cluster, histological subtype, TNM stages, clinicopathological stages, histological grade, and TP53 mutation, KRAS mutation, BRAF mutation and EGFR mutation as the annotations.

Multivariate survival analysis

To describe how these RNAs affected the prognosis of the BLCA patients, multivariate survival analysis was performed on the 511 signatures composed of the selected 9 biomarkers using permutation and combination method in the training dataset by using the “survival” package in the R language. Subsequently, we established a weighted overall survival (OS) prognostic index algorithm model for prediction of the prognosis using the following:

Risk score (RS) = ${\sum }_{i=1}^N\left(Exp\mathrm{\ast }Coef\right)$

where N was the number of prognostic PCGs, lncRNAs and microRNAs in the model, Exp stands for the expression value of the PCGs, lncRNAs and microRNAs, and Coef was the estimated regression coefficient of the PCGs, lncRNAs and microRNAs in the multivariate Cox regression model. Patients who have higher risk scores are expected to have a higher probability of a poor outcome. Selecting the median risk score in each dataset as a cutoff value, bladder urothelial carcinoma patients were divided into high- and low- risk groups (Zhou et al., 2015). Meanwhile, we performed Kaplan–Meier survival analyses to inspect the differential for survival distributions in different groups for each BLCA cohort, and the two-side log-rank test was used to assess the statistical significance. Furthermore, multivariate Cox regression analysis was conducted to test whether the risk score was independent of other clinical covariates.

Functional annotation via GO, KEGG analyses and gene set variation analysis (GSVA)

To further investigate the underlying biological roles and pathways of the three-dimension transcriptome signature, the co-expression relationships of the three PCGs, one lncRNA and one miRNA with the corresponding co-expressed protein-coding genes were calculated using Pearson correlation coefficients in the training dataset. GO analyses were conducted by using R package clusterProfiler (Yu et al., 2012). GSVA (Hanzelmann, Castelo & Guinney, 2013) was applied to obtain the abovementioned co-expressed PCGs based on the following criterion: Pearson correlation coefficient > —0.3—, p-value < 0.05.

Statistical analysis

Survival analysis was performed using the Kaplan–Meier method, and the differences between survival curves was assessed using the log-rank test. Univariate and multivariate analyses were conducted using Cox proportional hazard models. Survival predictive accuracy of prognostic models was estimated based on a receiver operating characteristic curve (ROC) analysis via “pROC” R package (Robin et al., 2011). All statistical analyses were performed in the R platform (R Core Team, 2019). All statistical tests were two-sided and p-values < 0.05 was considered statistically significant.

BLCA patients’ characteristics in the training and test sets

A total of four hundred patients included in this study were pathologically diagnosed with bladder urothelial carcinoma after preprocessing. Patients with missing overall survival data were excluded from this study. These patients were divided into a training set (n = 200) and a test set (n = 200) randomly and evenly. The patient age was classified as ≤ 61 and according to the X-tile software (Yale University Version 3.6.1, http://tissuearray.org) (Camp, Dolled-Filhart & Rimm, 2004). All of the baseline demographic and clinical characteristics of these three datasets are summarized in Table 1.

Identification of significant PCGs, lncRNAs and microRNAs associated with overall survival from the training set

Through preprocessing and preliminary screening of the significantly changed PCGs, lncRNAs and miRNAs in the training set (see methods), 8969 PCGs, 89 lncRNAs and 582 microRNAs expression profiles were generated for key RNAs identification.

We utilized three algorithms to identify the central RNAs used to construct the candidate prognostic prediction models. First, we conducted a univariate Cox proportional hazards regression analysis with the PCGs, lncRNAs and microRNAs expression profiling data as the independent variables while survival time and survival status were the dependent variables, and we identified 785 PCGs, 5 lncRNAs and 34 microRNAs that were significantly associated with patients’ OS (Table S1). Second, we used the robust likelihood-based survival models to identify 19 PCGs, 4 lncRNAs and 3 microRNAs that were the most highly correlated with the prognostic information. Third, using the random forests for survival regression and classification algorithm, we finally screened out 3 PCGs, 3 lncRNAs and 3 microRNAs according to the variable importance (VIMP) (Fig. S1).

Acquisition and construction of the prognostic PCG-lncRNA- microRNA signature model in the training set

The 3 PCGs, 3 lncRNAs and 3 microRNAs mentioned above in the training set could have 511 combinations and the corresponding risk score could be calculated according to every risk score model (Table S2). All of the risk scores of every patient were computed in the “Risk score (RS)” formula described in the methods. Simultaneously, we conducted ROC analysis for selecting a better prognostic signature. Finally, the PCG-lncRNA-microRNA model composed of ANXA1 (PCG), TPST1 (PCG), PSMB10 (PCG), DLEU1 (lncRNA) and miR-497-5p (miRNA) with a higher AUC and the significant regression coefficients (p-value < 0.05) was retained as the ultimate model (Table 2). The risk score of each patient was obtained with the criteria outlined below: Risk score = (0. 3024 × expression value of ANXA1) + (0. 1494 × expression value of TPST1) + (−0. 4377 × expression value of PSMB10) + (0. 4875 × expression value of DLEU1) + (0. 4262 × expression value of miR-497-5p). The AUC value of the PCG-lncRNA-microRNA signature mentioned above was 0.7669 (Fig. 2A), indicating that the model consisting of the five biomarkers could have a good performance for survival prediction. We also compared the survival predictive ability of the PCG-lncRNA-microRNA signature with pathologic stage using ROC analysis in the training dataset and concluded that the predictive power of the signature constructed was better than stage (AUC_signature = 0.7669, AUC_stage = 0.6591; p-value < 0.05) (Fig. 2B). In this study, we also compared the proposed model with other known prognostic signatures, including a six-gene signature (Wang et al., 2019) and an eight-mRNA signature (Zhu et al., 2019). As shown in the Figs. 2C–2D, the AUC value of the model proposed in this study was larger than the studies mentioned above. The result indicated that the prognostic signature model was candidate for predicting patients’ overall survival status. The results further indicated that the model constructed in our study was a novel predictive prognostic signature with high sensitivity and specificity in its clinical significance. Furthermore, from the “Risk score” formula, it suggested that ANXA1, TPST1, DLEU1 and miR-497-5p were possible risk factors and PSMB10 was a possible protective factor for survival (Table 2).

The 200 BLCA patients in the training set were assigned to the high-risk group (n = 100) and low-risk group (n = 100) according to the median risk score. Kaplan–Meier analysis indicated that patients in the high-risk group had a significantly poor outcome than those in the low-risk group (log-rank test, p-value < 0.0001; Fig. 3A). The distribution of risk scores, OS, vital status, and corresponding RNAs expression profiles of the 200 patients in the high-risk and low-risk groups are shown in Fig. 3B.

Validation of the prognostic signature model in the test set

Using the same method as in the training set, the patients in the test set (n = 200) and entire set (n = 400) were also classified into high-risk and low-risk groups to validate the survival prediction of the PCG-lncRNA-microRNA signature. In the test set, the 200 patients were divided into the high-risk group (n = 100) and low-risk group (n = 100). Kaplan–Meier curve analysis showed that patients in the high-risk group had a shorter OS than those in the low-risk group (log-rank test, p-value < 0.05; Fig. 3C). The AUC value of the signature was 0.645 (Fig. 2E). In the entire set, similar results were found (log-rank test, p-value < 0.0001; Fig. 3E). The AUC of the signature was 0.710 (Fig. 2F). However, the entire set was just used for exploration purposes instead of evaluating the value for the model. The distribution of risk scores, OS, vital status, and corresponding RNAs expression profiles of patients in the test set and entire set in the high-risk and low-risk groups are shown in Fig. 3D and Fig. 3F. We also noticed that ANXA1, TPST1, DLEU1 and miR-497-5p showed a tendency towards high expression in patients in the high-risk group, and PSMB10 presented with low expression in the high-risk group. These results are in accordance with the results in the training set (please see the section above) (Figs. 3D–3F).

The survival predictive ability of the three-dimension transcriptome signature is independent of other clinical features

To evaluate whether the PCG-lncRNA-microRNA signature maintained its prognostic power in the context of other clinical features, the results of multivariate Cox regression analysis showed that the power of the three-dimension transcriptome signature was maintained in survival predication and was significantly independent of other clinical features in the training set (hazard ratio [HR] = 4.209, 95% confidence interval [CI] [2.4541–7.218], p-value < 0.0001), the test set (HR = 1.679, 95% CI [1.0679–2.641], p-value < 0.05), and the entire set (HR = 2.4419, 95% CI [1.7262–3.454], p-value < 0.0001; Fig. 4).

NMF and SNF algorithm identifies three subtypes in BLCA

NMF was applied to the training set and cophenetic correlation coefficients were calculated to choose the appropriate number of clusters. Ultimately, a factorization rank of 3 was determined for the clusters by the method mentioned by Gaujoux & Seoighe (2010) who suggested choosing the smallest value of r for which this coefficient starts decreasing (Fig. S2A). The consensus matrix heatmap showed the preferable sharp boundaries, which indicated robust and stable clustering for the samples (Fig. S2B).

Then, 172 PCGs, 15 lncRNAs and 42 miRNAs metagenes identified by NMF were described as features, together with the 200 samples in the training set to build a similarity network to cluster cancer subtypes based on three-dimension transcriptome data. We classified all of the 200 BLCA patients into three clusters: cluster_1 (63 patients, 31.5%), clusters_2 (85 patients, 42.5%) and clusters_3 (52 patients, 26.0%). Statistical significance analysis of clustering showed that cluster_1/cluster_2, cluster_1/cluster_3 and cluster_2/cluster_3 had significant differences (p-value = 0) (Fig. S2C). The same classification method was performed in the test set and the entire set. Patients in these two datasets can also be divided into three clusters that presented with similar proportions (Fig. S2D). Furthermore, we found that patients in the high-risk group with high grade and stage I; and stage II; were enriched in subtype 2 and subtype 3 and few were in subtype 1, while patients in the low-risk group were enriched in subtype 1 and subtype 2 and few were in subtype 3. This suggested that subtype 2 and subtype 3 have potentially higher carcinogenic biological processes (Figs. S2E–S2G).

The average silhouette width between these clusters was 0.47 (range, from 0.42 to 0.55), which indicated the robustness of the classification in the training set (Fig. 5A). However, the overall survival probability between the three clusters had no significant differences (p-value = 0.2) (Fig. 5B). Meanwhile, the high-risk group had a poor prognosis compared with the low-risk group in every cluster (cluster_1: p-value < 0.0001, cluster_2: p-value < 0.0001, cluster_3: p-value = 0.00049) (Fig. S3). In the test set, the average silhouette width between the three clusters was 0.39 (range, from 0.25 to 0.54), the OS probability across the three clusters had no significant differences (p-value = 0.13) (Figs. 5C–5D), and the survival probability between the high-risk group and low-risk group had no significance except for the cluster_2 (p-value = 0.0085) (Fig. S3). In the entire set, the average silhouette width across the three clusters was 0.44 (range, from 0.31 to 0.57), the entire probability across the three clusters had significant differences (p-value = 0.015) (Figs. 5E–5F), while the high-risk group had a poor prognosis compared with the low-risk group in every cluster (cluster_1: p-value < 0.0001, cluster_2: p-value = 0.00059, cluster_3: p-value < 0.0001) (Fig. S3). The heatmap was presented in Fig. S4. The basal squamous, neuronal, luminal papillary and luminal infiltrated subtypes shared in the same cluster 3, which showed higher TP53 mutation and more dead events. And the majority of basal squamous and neuronal subtypes were enriched in the cluster 3, indicating high risk group, while the luminal papillary subtype was enriched in the cluster 2. Furthermore, the lncRNA and microRNA clusters also presented preferable general boundaries.

We also investigated the basic expression of these five biomarkers in the tumor and normal tissues (Fig. S5A). Based on the biomarkers of the constructed signature, we also investigated the potential biological significance behind these molecules. We calculated the expression correlation of the 5 genes, and the correlation coefficients of most of the genes were low in the subtypes or tumour samples. The genes with positive regression coefficients had positive correlations with each other. The relationship of ANXA1 and TPST1 showed a positive correlation (r = 0.21, p-value < 0.05) and the association of ANXA1 and DLEU1 also showed a positive correlation (r = 0.21, p-value < 0.05) (Fig. S5B). All of these indicated that these genes carried less overlapped information and showed low redundancy.

Functional characterization of GO and KEGG analyses and gene set variation analysis (GSVA)

To explore the co-expression relationships of the selected PCGs, lncRNAs and miRNAs with the PCGs, Pearson correlation coefficients were computed according to the standard. The expression of 4469/8969 protein-coding genes in the training set were highly correlated with that of at least one of the biomarkers. GO analysis, including biological process, cellular component and molecular function, was conducted. Focal adhesion and MAPK signaling pathway were also identified in the KEGG analysis. Then, GSVA for these co-expressed PCGs was performed based on the Hallmark gene sets from the Molecular Signatures Database (MSigDB). Several functionally related terms were identified and suggested that the selected five signatures might be involved in tumourigenesis by interacting with related PCGs that referred to important biological processes such as “epithelial to mesenchymal transition” (EMT), “KRAS signaling up” and so on (Fig. 6; Table S3).

Furthermore, the protein of the prognostic markers shared in the same subcellular localization, including plasma membrane, extracellular, nucleus and endosome and so on (Table S3). Data were from a subcellular localization database COMPARTMENTS (https://compartments.jensenlab.org/Search).