The RNA sequencing data (HTSeq-Counts), MuTect2-based somatic mutation data, and the corresponding clinical data of 412 urinary bladder cancer patients were downloaded from The Cancer Genome Atlas (TCGA) website (https://portal.gdc.cancer.gov/). The selection criteria of these patients included; (1) MIBC tissue samples; (2) having TP53 mutation data; (3) having gene expression data; (4) having follow-up data and survival status. A total of 402 MIBC patients were enrolled in the TCGA cohort. Then, HTSeq-Counts data were converted to transcripts per kilobase million (TPM) data and were log2-transformed (log2TPM) for subsequent analysis, which showed more comparability between samples (24). This study used the average expression value of genes in instances where gene duplication was detected.

The matrix files of the gene expression profile and clinical information of GSE32894, GSE48075, and GSE52219 data sets were extracted from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). Gene expression data of GSE48075 was log2-transformed. Notably, the average gene expression was used if multiple probes matched a single gene. Batch effects were removed through the “sva” package (version: 3.36.0; http://bioconductor.org/packages/release/bioc/html/sva.html). Then, GSE32894 (MIBC patients: 51) and GSE48075 (MIBC patients: 73) based on the same platform GPL6947 with follow-up data and survival status were integrated into the GEO cohort (n = 124). Twenty-three MIBC patients with information regarding their response to neoadjuvant chemotherapy (NAC) (NAC type: MVAC=methotrexate, vinblastine, adriamycin, and cisplatin) in the GSE52219 data set were included into this study. Data from the present study are publicly available from The Cancer Genome Atlas (TCGA) and GEO databases.

To identify potential differences in immunological pathways between TP53^mut (n=196) and TP53^wt (n=212) MIBC patients, GSEA (Version: 4.0; http://software.broadinstitute.org/gsea/index.jsp) was performed in TCGA cohort based on the reference gene set “c5.bp.v7.1.symbols.gmt.” The threshold was set at |Normalized Enrichment Score (NES)| > 1 and P < 0.01.

TMB is defined as the total number of deletions, insertions, base substitutions, or somatic gene coding errors detected per million bases. To explore the association between TP53 mutation and TMB, Perl scripts (version 5.30.2; http://strawberryperl.com/) were used to calculate the mutation rate with number of variants/the length of exons (38 million) for each sample of the TCGA cohort based on MuTect2-based somatic mutation data.

To identify DEIGs between TP53^mut (n=196) and TP53^wt (n=212) MIBC patients in the TCGA cohort, the “edgeR” package (version 3.30.0; http://www.bioconductor.org/packages/release/bioc/html/edgeR.html) was used with thresholds being false discovery rate (FDR) < 0.05 and |log2-fold change (FC)| > 1.0. The list of immune-related genes was downloaded from the ImmPort database (http://www.immport.org), which contains 17 immune categories based on different molecular function, such as TNF family receptors, B cell receptor signaling pathway, T cell receptor signaling pathway, and cytokines (25).

Gene Ontology (GO) enrichment analysis, including biological processes (BP), cellular components (CC), and molecular functions (MF), and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were conducted using “clusterProfiler” package (version 3.16.1; http://www.bioconductor.org/packages/release/bioc/html/clusterProfiler.html) (26).

Univariate Cox regression analysis was performed using “survival” package (version 3.1-12; https://cran.r-project.org/web/packages/survival/index.html) to assess the relationship between DEIGs and OS of MIBC patients in the TCGA cohort. DEIGs with P < 0.05 were considered prognostic immune-related genes. Further, the least absolute shrinkage and selection operator (LASSO) analysis was performed using the “glmnet” package (Version: 4.0; https://cran.r-project.org/web/packages/glmnet/index.html) to further screen prognostic DEIGs, which can resolve the collinearity problem and overfitting problem. Finally, TIPS was constructed by a multivariable Cox regression model based on the result from LASSO analysis, and the risk score of each patient was calculated by weighted estimators corresponding to the expression level of each gene. The patients in the TCGA cohort were divided into low- and high-risk prognostic groups as per the best cut off value determined by X-tile 3.6.1 software. To obtain the same formula and the uniform cutoff value to divide MIBC patients into low- and high-risk prognostic groups, a normalization for expression values of TIPS genes was conducted in the TCGA cohort and the GEO cohort with standard deviation (SD) = 1 and mean value = 0. To validate the TIPS, the risk score for each patient of the GEO cohort was calculated using the same formula and the patients were stratified into low- and high-risk prognostic groups following the same cutoff value obtained from the TCGA cohort. Of note, Kaplan-Meier survival analysis with the log-rank test was used to evaluate the survival difference between the low- and high-risk prognostic groups. Area under the curve (AUC) of time-dependent receiver operating characteristic (tROC) curve and the discrimination was used to evaluate prognostic performance of TIPS. The discrimination of TIPS was determined by a concordance index (C-index) using 1000 bootstrap resamples.

GSEA was performed to explore the enriched KEGG pathways of the TIPS genes. According to the median expression of each gene, 408 MIBC samples were divided into low- and high-expression groups. The gene set “c2.cp.kegg.v7.1.symbols.gmt” was chosen as the reference set, which was obtained from the Molecular Signatures Database V7.1 (MSigDB). |NES| > 1 and P < 0.01 were considered as statistically significant.

To investigate the association between TIPS and TILs, single sample gene set enrichment analysis (ssGSEA) was performed to determine the relative abundance of 28 subpopulations of TILs in the MIBC immune infiltrates, including Tregs, MDSC, TAMs, CD56^bright NK cells (NK) cells, etc (27, 28).

Notably, the immune checkpoints are the biomarkers for selecting patients with MIBC for immunotherapy. Therefore, this study analyzed the correlation between TIPS and critical immune checkpoints (PD1, CTLA4, LAG3, HAVCR2, and TIGIT).

Chemotherapy is effective for treating MIBC. Based on the Genomics of Drug Sensitivity in Cancer (GDSC) website, the clinical response to chemotherapy of each MIBC patient was estimated to explore whether there were differences in the response of low- and high-risk prognostic groups to chemotherapy and immunotherapy. Two commonly used chemotherapy drugs, gemcitabine and cisplatin, were selected for the chemotherapeutic response prediction, which involved the ridge regression based on the “pRRophetic” R package (https://github.com/paulgeeleher/pRRophetic20) to predict the half-maximal inhibitory concentration (IC50) of each TCGA-MIBC patient and the 10-fold cross-validation based on the GDSC training set was used to assess the prediction accuracy (29). To validate the feasibility of the TIPS to predict response to NAC, 23 MIBC patients with NAC in GSE52219 data set were included into this study. A normalization for the expression of TIPS genes was conducted in GSE52219 data set with a standard deviation (SD) = 1 and mean value = 0. The risk score for each patient was calculated using the same formula.

The immune checkpoint inhibitor has evolved to be the most potent tool for cancer therapy, such as anti-PD1 and anti-CTLA4, which activate the immune system to play an anti-tumor role (30). Here, subclass mapping on GenePattern website (https://cloud.genepattern.org/gp) was used to predict the anti-PD1 and anti-CTLA4 response of each MIBC patient as described previously (31, 32).

Out of the 402 MIBC patients in the TCGA cohort with follow-up and survival status data, 342 MIBC patients with complete clinical data, including age, gender, T stage, N status and TMB, were selected for subsequent analyses. Univariate and multivariate Cox regression analysis was performed to identify whether the prediction of TIPS was independent from other clinical characteristics (age, gender, T stage, N status and TMB). To further identify whether the TIPS was a statistically significant prognostic factor regardless of clinical data, the patients were stratified based on age, T stage and TMB, and Kaplan-Meier survival analysis was used to assess the difference in the OS between low- and high-risk prognostic groups of different subgroups.

To facilitate the prediction of 1-, 3-, and 5-year OS probability in MIBC patients, a nomogram was developed based on the results from multivariate Cox regression analyses using the “rms” R package (version 6.0-0; https://cran.r-project.org/web/packages/rms/index.html). The commonly used methods i.e., AUC of tROC, C-index using 1000 bootstrap resamples, and calibration plot were used to validate the performance of the nomogram. Calibration plots were used to assess the consistency between nomogram-predicted probabilities and observed probabilities using 500 bootstrap resamples, with the 45degree line representing the ideal predicted values. Moreover, decision curve analysis (DCA) was used to determine the net benefits derived from the use of the nomogram, TIPS, age, T stage, N status and TMB (33).

The 95% confidence interval (CI) and HR (Hazard Ratio) were generated by Cox regression analysis and Kaplan-Meier survival analysis. Cox regression analysis assume that the HR is constant over time, therefore, is significant to assess the validity of the proportional hazards assumption (PH assumption). The scale Schoenfelder residual test based on the “survival” R package (Version 3.2-7; https://cran.r-project.org/web/packages/survival/index.html) was performed to evaluate whether clinical variables violated the PH assumption (34). The FDR method was used to control for multiple testing, which is the percentage of important tests that can lead to a false positive (35). Spearman correlation analysis was performed to access the existence of a correlation between variables. Mann-Whitney-Wilcoxon Test or Student’s t test was used to evaluate the comparison between groups for continuous variables. Categorical variables were compared between groups using Chi-square or Fisher’s exact tests. All statistical analyses were performed in R software (version 3.6.2; https://www.r-project.org/) and they were 2-sided, with a P-value of less than 0.05 considered statistically significant.

As shown in Figure 1A , TP53 mutation was the most frequent type in TCGA-MIBC cohort (47%), with missense mutation being the main type. In order to determine the prognosis of different types of TP53 mutation, we conducted stratification analysis based TP53 patterns which combined the nonsense mutation, frameshifts (deletions and insertions), in frame shift deletion, splice site and multi-hit mutations to non-missense mutations. The Kaplan-Meier survival analysis suggested that patients with missense mutations have a poorer prognosis when compared with non-missense mutations ( Figure 1B ). Previous studies found that TP53 mutation is associated with the pathological staging and prognosis of MIBC. However, research on the association between TP53 mutation and immunophenotype in MIBC has not matured. Therefore, for the first time, this study used RNA sequencing data and clinical information of TCGA-MIBC patients to identify immune-related biological processes related to TP53 status. GSEA analysis showed that TP53^mut MIBC patients were significantly enriched in 657 biological processes, out of these, five were immune-related biological processes ( Figure 1C ). In contrast, no immune-related biological processes were enriched in the TP53^wt MIBC patients. In addition, the proportion of TP53^mut MIBC patients was significantly high among high-TMB group compared with low-TMB group (Chi-square test, P value <0.001; Figure S1A ). The TMB of TP53^mut MIBC patients was higher than TP53^wt TMB patients (Student’s t test, P value = 0.038; Figure S1B ), suggesting that the TP53 mutation was closely associated with high TMB.

Based on the results obtained from the GSEA analysis, the TP53 mutation was closely related to immune-related biological processes, hence, MIBC patients were subdivided into TP53^mut and TP53^wt groups, then DEIGs were explored to further identify the correlations between the TP53 mutation and immunophenotype in MIBC. A total of 50 upregulated genes and 49 downregulated genes were identified (FDR < 0.05 and |log2- FC| > 1.0) ( Figures 2A, B ). The differentially expressed gene analysis using the “edgeR” package is shown in supplement Table S1 . As shown in Figure S2A , significantly enriched BP of DEIGs were detected, including B cell proliferation, second messenger mediated signaling and humoral immune response. Several CC GO terms were detected, including external side of plasma membrane, secretory granule lumen and cytoplasmic vesicle lumen ( Figure S2B ). In GO terms of MF, receptor ligand activity, signaling receptor activator activity and cytokine activity were significantly enriched terms ( Figure S2C ). According to the KEGG pathway analysis, cytokine-cytokine receptor interaction, IL−17 signaling pathway and natural killer cell mediated cytotoxicity were mostly associated with the DEIGs ( Figure S2D ).

Univariate Cox regression analysis was performed, where 13 of the 99 DEIGs were identified to be significantly associated with OS in TCGA-MIBC patients ( Table S2 ). For further screen the 13 prognostic DEIGs, LASSO analysis was applied ( Figures 3A, B ). The LASSO analysis was used to shrink all regression coefficients toward zero and to select variables simultaneously; and the optimal lambda value were determined through 10-fold cross-validations. Finally, we identified TIPS using multivariate Cox regression analysis that was significantly associated with OS in MIBC patients. The risk score of each TCGA-MIBC patient was calculated as follows: [−0.42 × Expression value of ORM1] + [0.17 × Expression value of PTHLH] + [−0.39 × Expression value of CTSE] ( Table S3 ). Figure 2B showed the differential expression of these three genes between TP53^mut and TP53^wt MIBC patients. Further, the TCGA-MIBC patients were divided into either low- or high-risk prognostic group based on −0.037 determined by X-tile 3.6.1 software ( Figures S3A, B ). The risk score distribution, gene expression, and OS status of TCGA-MIBC patients were shown and ranked based on risk score values of the TIPS ( Figure 3C ).

In addition, Kaplan-Meier survival analysis showed that high-risk prognostic group with the higher risk score had a significantly poorer prognosis than low-risk prognostic group with the low-risk score (HR = 2.22, 95% CI = 1.62–3.06, P < 0.0001) ( Figure 3D ). The AUCs for 1-, 3-, and 5-year OS predictions for TIPS were 0.625, 0.643, and 0.640 respectively, indicating that TIPS had a satisfactory sensitivity and specificity ( Figure 3E ). C-index of the TIPS was 0.62 (95%CI: 0.57–0.66, P<0.0001). As shown in Figures S4A to C , the AUCs for 1-, 3- and 5-year OS predictions of TIPS were higher than AUCs of age, gender, T stage, N status and TMB in the TCGA cohort.

To evaluate the robustness of TIPS, its performance was assessed in the GEO cohort, which consisted of 124 MIBC patients. The patients in the GEO cohort were stratified into low- and high-risk prognostic groups using the same formula and the same cutoff value (−0.037) obtained from the TCGA cohort ( Table S4 ). The risk score distribution, gene expression, and OS status of patients in the GEO cohort are shown in Figure 3F . In comparison with the high‐risk prognostic group in the GEO cohort, a significantly higher survival rates were observed in the low-risk prognostic group (HR = 2.05, 95% CI= 1.26–3.33, P = 0.0084) ( Figure 3G ). This was in line with the results of the TCGA cohort. Moreover, the TIPS achieved an AUC of 0.650 at 1 year, 0.691 at 3 years, and 0.671 at 5 years ( Figure 3H ). C-index of the TIPS was 0.61 (95%CI: 0.55–0.67, P<0.0001). The above results demonstrate the TIPS is effective in predicting the OS of MIBC patients.

To further investigate enriched KEGG pathways of CTSE, ORM1 and PTHLH in MIBC, GSEA was performed based on TCGA-MIBC RNA-seq data. As shown in Figure S5 , genes in the high expression groups of CTSE (A) were mainly enriched in the metabolism-related pathway, such as “alpha linolenic acid metabolism,” “arachidonic acid metabolism” and “drug metabolism cytochrome P450.” Meanwhile, the “Natural killer cell mediated cytotoxicity,” “complement and coagulation cascades,” and “cytokine-cytokine receptor interaction” were enriched in the high expression groups of ORM1 (B) and PTHLH (C), whereas “T cell receptor signaling pathway” and “Nod like receptor signaling pathway” were enriched in ORM1 and PTHLH high-expression groups, respectively. The results of GSEA revealed that ORM1 and PTHLH correlated with immune signaling pathways.

Using the ssGSEA method, the differences in 28 subpopulations of TILs in the MIBC immune infiltrates between low- and high-risk prognostic groups were identified. Figures 4A, B are a summary of the results obtained from the 402 MIBC patients, showing that multiple subpopulations of TILs were significantly different between low- and high-risk prognostic groups. Among them, it was worth noting that the high-risk prognostic group had a significantly higher abundance of Tregs, TAM, and MDSCs as well as significantly lower abundance of CD56^bright NK cells compared to the low-risk prognostic group (P < 0.05). The heterogeneity of TILs in the MIBC immune infiltrates indicated that the TILs may be served potential prognostic biomarkers, immunotherapy targets, and might exhibit important clinical significance. Besides, it was found that the risk score was positively related to the expression of critical immune checkpoints in the TCGA cohort, including PD1, CTLA4, LAG3, HAVCR2, and TIGIT (P<0.05) ( Figure 5A ). Also, the differential expression of PD1, CTLA4, LAG3, HAVCR2, and TIGIT between the low- and high-risk prognostic groups were analyzed. Results showed that the expression of PD1, CTLA4, LAG3, HAVCR2, and TIGIT was significantly higher in high-risk prognostic group than in low-risk prognostic group ( Figures 5B–F ), implying that the bad prognosis of high-risk prognostic group might be associated with the immunosuppressive microenvironment.

The high APOBEC-signature mutation load may stimulate a natural host immune reaction to curb tumor growth and metastasis and was associated with better prognosis of MIBC (36). The expression of APOBEC3A and APOBEC3B positively correlated with levels of APOBEC-signature mutagenesis, which was ubiquitous and carcinogenic (36, 37). The file including the information of the APOBEC-signature mutation load was downloaded from the study of Robertson et al. (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5687509/#) to analyze the association between TIPS and APOBEC-signature mutation. Low-risk prognostic group had a high APOBEC-signature mutation load and low expression of APOBEC3A and APOBEC3B, but this was not statistically significant ( Figures S6A–C ).

Given that chemotherapy is a common approach in the treatment of MIBC, this study attempted to evaluate the response of low- and high-risk prognostic groups to the treatment by gemcitabine and cisplatin. The IC50 of each TCGA-MIBC patient was estimated and as shown in Figures 6A, B , it resulted in a significant difference in the estimated IC50 between low- and high-risk prognostic groups for the two chemotherapy drugs, where the high-risk prognostic group was more sensitive to the treatment (P = 0.009 for cisplatin, P = 0.001 for gemcitabine). As shown in Figure 6C , the proportion of NAC responders was significantly high among high-risk prognostic group compared with the low-risk prognostic group in GSE52219 data set (Fisher’s exact test, p value < 0.001). High-risk prognostic group were more sensitive to NAC and may benefit from it which was consistent with the above predictions. To determine the potential of the risk score as a biomarker for immunotherapy, the difference in the estimated immunotherapeutic response of low- and high-risk prognostic groups was identified. Interestingly, the analysis suggested that high-risk prognostic group might be more responsive to anti–PD-1 therapy (Bonferroni corrected P = 0.023) ( Figure 6D ).

Baseline characteristics of 342 MIBC patients in the TCGA cohort are summarized in Table 1 . TMB was stratified into low and high TMB groups according to 5.55/Mb determined by X-tile 3.6.1 software ( Figure S7 ). Regarding the Schoenfeld Individual Test, the P value of five clinical characters, including age (A), gender (B), T, stage (C), N status (D), TIPS (E) and TMB (F) was >0.05 and the PH assumption could not be considered violated ( Figure S8 ). Univariate and multivariate Cox regression analyses were performed to explore the independence of TIPS to clinical characteristics in the TCGA cohort. After adjustment for clinical characteristics and multiple testing using FDR method, the TIPS remained an independent prognostic signature indicating its robustness in independently predicting the OS of MIBC patients ( Figures 7A, B ). Subsequently, Kaplan-Meier survival analysis was used to analyze the predictive value of the TIPS in different subgroups stratified by age, T stage and TMB, showing that the high-risk prognostic group had a poorer OS compared to low-risk prognostic group in ≤65 years, >65 years, T2, T3&T4, low and high TMB subgroups ( Figures 7C–H ). We further analyzed the relationship between the risk score and the clinical characteristics of MIBC (including age, T stage and N status). In terms of T stage, T3&T4 stage patients had higher risk scores than T2 stage patients in the TCGA data set but the risk scores for the N1 status and >65 years patients were not higher than those for the N0 status and ≤65 years patients ( Figures S9A–C ).

To facilitate the clinical decision making, a nomogram was developed based on the TIPS to predict 1-, 3-, and 5-year OS of the TCGA-MIBC patients. The risk score, age, T stage, N status and TMB were selected into the nomogram by a stepwise Cox regression model ( Figure 8A ). And as shown in Figure S4A-C , the AUC for 1-, 3- and 5-year OS predictions of the nomogram was higher than AUC of the TIPS, age, gender T stage, N status and TMB in the TCGA cohort. The C-index of the nomogram was 0.69 (95% CI, 0.64–0.73), while that of the risk score was 0.62 (95% CI, 0.57–0.66) in the TCGA cohort and 0.61 (95% CI, 0.55–0.67) in the GEO cohort. The calibration plots showed consistency in predicting OS of the nomogram with the actual probability of OS at 1-, 3- and 5-year ( Figure 8B ). Moreover, the DCA for 5-year OS prediction showed that the nomogram had the highest net benefit across 0% to 80% threshold probabilities ( Figure 8C ). Meanwhile, the net benefit of the TIPS was higher than age, T stage, N status and TMB.