Development of stemness-related gene signature
Firstly, a total of 25 prognostic genes were identified by univariate survival analysis (Additional file 2: Table S2). Among these genes, six genes were regarded as favorable genes for OS (HR < 1.00) while other 19 genes were unfavorable genes (HR > 1.00). With LASSO regression (Fig. 1a/b) and stepwise multivariate Cox regression, an optimal seven-gene risk signature (Table 2) was constructed, which included five unfavorable genes (EGFR, FOXA2, MME, SMOC2, and TFRC) and two favorable genes (HES1, RBM6). The risk score of each patient was generated as follow: risk score = (EGFR*0.988) + (FOXA2*1.721) + (HES1*-1.817) + (MME*0.797) + (RBM6*-3.190) + (SMOC 2*0.916) + (TFRC*1.711). In training cohort, 333 patients were classified into high-risk group (n = 167) and low-risk group (n = 166) by the median risk score 1.051. Kaplan-Meier curve exhibited a significant difference in OS between the two groups as high-risk group of patients had poorer clinical outcome than low-risk group (p < 0.001, Fig. 1c). The time-dependent ROC curve suggested that the AUCs of the stemness-related signature were 0.673, 0.732, 0.760 at 1-year, 3-years, and 5-years respectively (Fig. 1d). In addition, the seven-gene signature had better sensitivity and specificity than clinical factors in predicting 1-year, 3-years, and 5-years OS in MIBC patients (Additional file 3: Figure S1).
Table 2 The stemness-related seven-gene risk signature.
Gene
|
Coefficient
|
HR
|
HR.95L
|
HR.95H
|
p
|
EGFR
|
0.988
|
4.677
|
1.778
|
12.302
|
0.002
|
FOXA2
|
1.721
|
2.873
|
1.013
|
8.146
|
0.047
|
HES1
|
-1.817
|
0.155
|
0.044
|
0.548
|
0.004
|
MME
|
0.797
|
2.993
|
1.383
|
6.477
|
0.005
|
RBM6
|
-3.190
|
0.192
|
0.047
|
0.780
|
0.021
|
SMOC2
|
0.916
|
2.280
|
1.081
|
4.809
|
0.030
|
TFRC
|
1.711
|
6.038
|
1.423
|
25.625
|
0.015
|
HR: hazard ratio. HR.95L/H: lower /upper limit of the 95% confidence interval of hazard ratio.
Independence of stemness-related gene signature
In univariate Cox regression analysis, the seven-gene risk signature was a strong risk factor for OS (Hazard ration (HR) = 1.919, p < 0.001, Table 3). Clinical factors including age (HR = 1.962, p < 0.001), overall stage (HR = 1.958, p < 0.001), T stage (HR = 1.652, p < 0.001), and N stage (HR = 1.589, p < 0.001) were also risk factors for OS (Table 3). Moreover, the seven-gene risk signature (HR = 1.925, p < 0.001, Table 3) was an independent prognostic indicator for OS in MIBC patients when adjusted for age, gender, and clinical stage. In training cohort, Kaplan-Meier survival analysis was conducted for patients stratified by age, gender, overall stage, T stage, and N stage. High-risk group of patients always suffered worse prognosis when compared with their low-risk counterparts in all clinical subgroups and the results were significant (Fig. 2). Taken together, the novel gene signature could independently predict OS in MIBC patients.
Table 3 Univariate and multivariate Cox regression analysis of the gene signature and clinical factors in training cohort.
|
Univariate analysis
|
|
|
Multivariate analysis
|
|
|
Variable
|
HR (95%CL)
|
p
|
|
HR (95%CL)
|
p
|
|
age
|
1.962 (1.326−2.901)
|
<0.001
|
|
1.969 (1.327−2.922)
|
<0.001
|
|
gender
|
0.871 (0.596−1.271)
|
0.473
|
|
0.922 (0.627−1.355)
|
0.677
|
|
stage
|
1.958 (1.526−2.513)
|
<0.001
|
|
1.582 (1.009−2.481)
|
0.046
|
|
T
|
1.652 (1.268−2.152)
|
<0.001
|
|
1.164 (0.839−1.615)
|
0.363
|
|
N
|
1.589 (1.330−1.899)
|
<0.001
|
|
1.140 (0.830−1.567)
|
0.418
|
|
Risk score
|
1.919 (1.593−2.312)
|
<0.001
|
|
1.925 (1.564−2.370)
|
<0.001
|
|
HR:hazard ratio, CL: confidence interval.
Validation of stemness-related gene signature
We further verified the performance of the seven-gene signature in three validation cohorts. Patients in GSE13507 (high risk = 69, low risk = 96) and GSE32548 (high risk = 62, low risk = 65) were categorized by the cut-off value 1.051 in training cohort. While in GSE48075, only 6 patients were low-risk and the remaining 66 patients were all high-risk according to the cut-off value 1.051. Therefore, we reassigned the 72 patients in GSE48075 into equal groups (n = 36, 36 respectively) by their own median risk score 5.850. Kaplan-Meier survival curves showed that high-risk patients also suffered reduced overall survival when compared with low-risk patients, which was observed in GSE13507 (p < 0.001, Fig. 3a), GSE32548 (p < 0.001, Fig. 3b) and GSE48075 (p < 0.001, Fig. 3c). In GSE13507, the AUCs of seven-gene signature were 0.732, 0.645, 0.629 at 1, 3, 5 years, respectively (Fig. 3d). In GSE32548, the AUCs of seven-gene signature were 0.773, 0.803, 0.709 at 1, 3, 5 years respectively (Fig. 3e). In GSE48075, the AUCs at 1, 3, 5 years were 0.729, 0.694, 0.608 respectively (Fig. 3f).
Association between stemness-related gene signature and clinical features
The relationship between the stemness-related gene signature and clinical factors was analyzed. As expected, the risk score was significantly elevated along with cancer progression, including advanced T stage (p < 0.001, Fig. 4a), N stage (p = 0.018, Fig. 4b), and overall stage (p < 0.001, Fig. 4c). However, no difference in risk score was observed between different subgroups of age and gender. Moreover, the seven-gene risk signature was closely linked to response to primary chemotherapy and radiotherapy as the risk score was sequentially increasing in patients with CR, PR, SD, and PD (p = 0.002, Fig. 4d). These findings suggested that high risk score based on seven-gene signature was generally associated with unfavorable clinical manifestations.
Comparation between stemness-related gene signature and mRNA stemness indices
The mRNAsi of TCGA bladder cancer patients were obtained from Malta et al.’s study. We matched the samples between clinical data and mRNAsi by sample id number. As a result, a total of 401 patients were provided with mRNAsi. We used the median mRNAsi to divide 401 patients into two groups. Kaplan–Meier survival curve showed that patients with higher mRNAsi had significantly poor long-term survival than patients with lower mRNAsi (p < 0.001, Fig. 4e). The time-dependent ROC curve demonstrated that the AUCs of mRNAsi were 0.559, 0.548, 0.597 at 1-year, 3-years, and 5-years respectively (Fig. 4f). Furthermore, we also identified a significant correlation between mRNAsi and our seven-gene signature. The high-risk group patients identified by our seven-gene signature had significantly higher mRNAsi (Fig. 4g).
Integrated nomogram by combining the stemness-related gene signature with clinical factors
We subsequently developed a nomogram in training cohort to achieve better predictive performance (Fig. 5a). Three independent risk factors for OS including age (HR = 1.969, p < 0.001), overall stage (HR = 1.582, p = 0.046), and the seven-gene risk signature (HR = 1.925, p < 0.001) were included in the nomogram (Table 3). No collinearity was observed between three risk factors. The AUCs of the nomogram were 0.761, 0.785, 0.776 at 1-, 3-, 5-years respectively (Fig. 5b). The calibration plot exhibited high degree of agreement between predicted and actual survival probabilities (Fig. 5c-e).
Gene Set Enrichment Analysis
Using GSEA, we explored the signaling pathways which were associated with our gene signature-based risk classification. Under the criteria of p < 0.05 and q < 0.05, a total of 1825 biological processes were significantly altered in high risk group. The most altered biological processes included extracellular structure organization (NES = 2.400), wound healing (NES = 2.390), cell matrix adhesion (NES = 2.367), ossification (NES = 2.340), angiogenesis (NES = 2.333), cartilage development (NES = 2.328), cell substrate adhesion (NES = 2.325), connective tissue development (NES = 2.318), vasculature development (NES = 2.318 ), and regulation of cell shape (NES = 2.317), which were dominantly related with extracellular adhesion and angiogenesis (Figure 6; Additional file 4: Table S3,).