Utility of a Molecular Signature for Predicting Recurrence and Progression in Non-Muscle-Invasive Bladder Cancer Patients: Comparison with the EORTC, CUETO and 2021 EAU Risk Groups

To evaluate the utility of different risk assessments in non-muscle-invasive bladder cancer (NMIBC) patients, a total of 178 NMIBC patients from Chungbuk National University Hospital (CBNUH) were enrolled, and the predictive value of the molecular signature-based subtype predictor (MSP888) and risk calculators based on clinicopathological factors (EORTC, CUETO and 2021 EAU risk scores) was compared. Of the 178 patients, 49 were newly analyzed by the RNA-sequencing, and their MSP888 subtype was evaluated. The ability of the EORTC, MSP888 and two molecular subtyping systems of bladder cancer (Lund and UROMOL subtypes) to predict progression of 460 NMIBC patients from the UROMOL project was assessed. Cox regression analyses showed that the MSP888 was an independent predictor of NMIBC progression in the CBNUH cohort (p = 0.043). Particularly in patients without an intravesical BCG immunotherapy, MSP888 significantly linked with risk of disease recurrence and progression (both p < 0.05). However, the EORTC, CUETO and 2021 EAU risk scores showed disappointing results with respect to estimating the NMIBC prognosis. In the UROMOL cohort, the MSP888, Lund and UROMOL subtypes demonstrated a similar capacity to predict NMIBC progression (all p < 0.05). Conclusively, the MSP888 is favorable for stratifying patients to facilitate optimal treatment.


Introduction
In 2022, an estimated 60,800 adults will be newly diagnosed with non-muscle invasive bladder cancer (NMIBC) in the United States, which is approximately 75% of all bladder cancer (BCa) patients, and in younger patients (<40 years of age) this percentage is even higher [1]. NMIBC is one of the costliest cancers to treat due to the heterogeneity of the disease and high recurrence rates; this means that patients must be monitored throughout their life [2]. Treatment of NMIBC aims mainly to reduce recurrence and to prevent progression to a muscle-invasive phenotype [3]. Therefore, predicting recurrence and progression in individual patients is crucial to facilitate optimal treatment and management protocols.

Construction of the MSP888 Classifier and Validation of Its Prognostic Value in the CBNUH RNA-seq Cohort
The molecular clusters in the Chungbuk National University Hospital (CBNUH) RNAsequencing (RNA-seq) cohort were defined based on the MSP888 classifier. As in our previous study, the optimal gene sets, comprising 888 genes, whose expression values were significantly diverse among Clusters 1, 2 and 3 were used to classify the molecular subgroups of NMIBC in the CBNUH RNA-seq cohort. A deep belief network-based (DBN) [14] deep learning method was then applied to develop and train the prediction model ( Figure 1A). Supplementary Table S1 shows the detailed data, along with predicted cluster results. Kaplan-Meier survival plots showing the recurrence-free survival (RFS) and progression-free survival (PFS) of NMIBC patients sub-grouped by the MSP888 prediction model were generated ( Figure 1B-Figure 2E). The recurrence and progression rates of NMIBC patients in Cluster 3 were higher than those in the other clusters ( Figure 1B,D); however, due to the small sample size with a low incidence rate, the p-Value was not significant. When we compared the RFS and PFS in individual clusters, the recurrence rate of NMIBC patients classified into Cluster 2 was dependent on BCG treatment that patients treated with BCG suffered worse RFS (p = 0.001) ( Figure 1C). of NMIBC patients classified into Cluster 2 was dependent on BCG treatment that patients treated with BCG suffered worse RFS (p = 0.001) ( Figure 1C). The ability of the MSP888 classifier to predict recurrence (C) or progression (D) after BCG treatment in three subgroups. Patients in Cluster 2 (E) derived significant benefit from BCG treatment. Data were plotted according to whether patients received BCG therapy. BCG, Bacillus Calmette−Guérin; CBNUH, Chungbuk National University Hospital; MSP888, a molecular signature−based subtype predictor comprising Cluster 1 (EP, equivocal prognosis), Cluster 2 (REC.BCG+, related to recurrence and response to BCG treatment), and Cluster 3 (DP.BCG+, related to progression and response to BCG treatment).

MSP888 Classifier Predicts NMIBC Prognosis in the CBNUH Cohort Better Than the EORTC, CUETO and 2021 EAU Risk Groups
In the previous study, we used 129 NMIBC samples to estimate the power of MSP888 for predicting NMIBC prognosis [13]. Here, we examined 49 cases in the CBNUH RNA−seq cohort and then compared the prognostic ability of the MSP888 classifier to that of EORTC, CUETO and 2021 EAU risk groups acquired from a total of 178 NMIBC patients. The correlation between each risk estimator and the prognosis of NMIBC patients was explored using Pearson's Chi−square test. The interval likelihood ratio (LR) was calculated as the probability of an individual test result occurring when recurrence/progression is present divided by the probability of an individual test result occurring when recurrence/progression is absent in order to demonstrate how likely it is that a patient will experience recurrence or progression

Analysis of NMIBC Recurrence
Only the MSP888 classifier showed a significant association with recurrence in patients without BCG treatment (Pearson's Chi−square p = 0.024) ( Table 1). When we calculated the interval LR for recurrence in the MSP888 classifier of NMIBC patients not treated with BCG, we found that Cluster 2 had a significantly higher LR of recurrence (LR, 1.942; 95% CI, 1.257-3.000) than Clusters 1 or 3 (Table 1). However, neither the EORTC nor the CUETO risk scores correlated with NMIBC recurrence. Kaplan-Meier survival plots demonstrated that patients in Cluster 2 suffered more recurrence than those in Cluster 1 (log−rank test, p = 0.027) ( Figure 2B). In addition, univariate and multivariate Cox regression analyses of the three clusters in the CBNUH cohort, along with known clinicopathological risk factors and the EORTC and CUETO risk scores, were conducted to determine Recurrence-free survival (B) and progression-free survival (D) of NMIBC patients are predicted by different clusters. The ability of the MSP888 classifier to predict recurrence (C) or progression (D) after BCG treatment in three subgroups. Patients in Cluster 2 (E) derived significant benefit from BCG treatment. Data were plotted according to whether patients received BCG therapy. BCG, Bacillus Calmette-Guérin; CBNUH, Chungbuk National University Hospital; MSP888, a molecular signature-based subtype predictor comprising Cluster 1 (EP, equivocal prognosis), Cluster 2 (REC.BCG+, related to recurrence and response to BCG treatment), and Cluster 3 (DP.BCG+, related to progression and response to BCG treatment).

MSP888 Classifier Predicts NMIBC Prognosis in the CBNUH Cohort Better Than the EORTC, CUETO and 2021 EAU Risk Groups
In the previous study, we used 129 NMIBC samples to estimate the power of MSP888 for predicting NMIBC prognosis [13]. Here, we examined 49 cases in the CBNUH RNAseq cohort and then compared the prognostic ability of the MSP888 classifier to that of EORTC, CUETO and 2021 EAU risk groups acquired from a total of 178 NMIBC patients. The correlation between each risk estimator and the prognosis of NMIBC patients was explored using Pearson's Chi-square test. The interval likelihood ratio (LR) was calculated as the probability of an individual test result occurring when recurrence/progression is present divided by the probability of an individual test result occurring when recurrence/progression is absent in order to demonstrate how likely it is that a patient will experience recurrence or progression

Analysis of NMIBC Recurrence
Only the MSP888 classifier showed a significant association with recurrence in patients without BCG treatment (Pearson's Chi-square p = 0.024) ( Table 1). When we calculated the interval LR for recurrence in the MSP888 classifier of NMIBC patients not treated with BCG, we found that Cluster 2 had a significantly higher LR of recurrence (LR, 1.942; 95% CI, 1.257-3.000) than Clusters 1 or 3 (Table 1). However, neither the EORTC nor the CUETO risk scores correlated with NMIBC recurrence. Kaplan-Meier survival plots demonstrated that patients in Cluster 2 suffered more recurrence than those in Cluster 1 (log-rank test, p = 0.027) ( Figure 2B). In addition, univariate and multivariate Cox regression analyses of the three clusters in the CBNUH cohort, along with known clinicopathological risk factors and the EORTC and CUETO risk scores, were conducted to determine prognostic factors. CUETO risk scores were calculated to determine prognostic factors. Univariate analysis identified tumor grade (both 1997 and 2004 WHO grading systems), the CUETO score and MSP888 as significant prognostic indicators of RFS in NMIBC patients not treated with BCG (all p < 0.05) ( Table 2). Because the tumor grade is an index used to calculate the CUETO score, we supposed that the prognostic significance of the CUETO score may be affected by tumor grade (1997 WHO grading system); therefore, we included only the tumor grade and the MSP888 classifier in multivariate analysis. Multivariate Cox analysis revealed that the MSP888 classifier retained statistical significance for RFS of NMIBC patients (HR, 4.104; 95% CI, 1.319-12.770; p = 0.015) ( Table 2). These results illustrate the high prognostic relevance of Cluster 2, characterized as a REC.BCG+ subtype associated with a poor prognosis for disease recurrence, but with a good response to BCG treatment.

Analysis of NMIBC Progression
Here, we found that the MSP888 classifier was able to differentiate NMIBC patients with progression from those without (Pearson's Chi-square, p < 0.0001). In addition, the interval LR suggested that patients in Cluster 3 were more likely to suffer progression than those in Clusters 1 or 2 (LR, 2.738; 95% CI, 1.863-4.024); however, none of the EORTC, CUETO and 2021 EAU risk scores correlated with NMIBC progression (Table 1). In particular, both the MSP888 classifier and the CUETO risk score showed a relevant correlation in NMIBC patients treated with BCG treatment: the LR for patients in Cluster 3 and for patients with a high CUETO risk score were 3.840 and 4.000, respectively, which is higher than for the other groups (Table 1). Kaplan-Meier survival plots revealed that the PFS of patients in Cluster 3 was poorer than that for patients in Cluster 1 or 2 (log-rank test, p < 0.0001, respectively) ( Figure 2D and E). Univariate Cox analysis identified age, tumor grade (both 1997 and 2004 WHO grading systems), the CUETO score and MSP888 as signif-icant indicators of PFS in NMIBC patients ( Table 3). As mentioned above, the prognostic significance of the CUETO score may be largely influenced by tumor grade (1997 WHO grading system); therefore, multivariate Cox analyses were evaluated using tumor grade and the MSP888 classifier. The results showed that the MSP888 classifier was a crucial independent risk factor for NMIBC progression (HR, 4.285; 95% CI, 1.708-10.750; p = 0.002); this was particularly true for patients not treated with BCG, in whom the HRs of the MSP888 classifier were 14.619 (95% CI, 4.042-52.869; p < 0.0001) and 7.411 (95% CI, 1.985-27.669; p = 0.003) ( Table 3). These results suggest the high prognostic value of Cluster 3, characterized as the DP.BCG+ subtype; this subtype has a poor prognosis with respect to disease progression, but a good response to BCG treatment. Taken together, the above results indicate that the MSP888 classifier is a more powerful predictor of NMIBC prognosis than the EORTC, CUETO and 2021 EAU risk scores.

In the UROMOL Cohort, MSP888 Classifier Predicts NMIBC Prognosis as well as Previously Developed Molecular Classifications
The UROMOL project is a European multicenter prospective study that aims to identify molecular markers that predict the likelihood of progression in patients with NMIBC [9]. The UROMOL study sub-grouped NMIBC into three major classes with basaland luminal-like characteristics and different clinical outcomes, which is partly consistent with Lund subtypes [9]. In the current study, we compared the predictive efficacy of the MSP888 classifier with that of the UROMOL classification, the Lund subtype and the EORTC score in 460 NMIBC patients enrolled in the UROMOL project. Due to a lack of clinical information recorded by the UROMOL project, we could only assess the risk of progression. Table 4 summarizes the associations between NMIBC progression and the four risk estimators: MSP888 classifier, UROMOL classification, Lund subtype and EORTC risk score. All four risk predictors correlated with NMIBC progression, especially in patients not treated with BCG (Pearson Chi-square p < 0.05). The interval LR revealed that class 2 of the UROMOL classification, the genomically unstable Lund subtype, a high EORTC risk score and Cluster 3 of the MSP888 were associated with a higher risk of progression. Thus, the predictive performance of the four risk estimators for NMIBC progression is in line with that reported previously in the UROMOL and CBNUH cohort studies. To identify and compare the prognostic independence of the MSP888 classifier, the UROMOL classification and the Lund subtype with respect to NMIBC progression, we applied Cox regression analyses to each risk predictor in the UROMOL cohort, along with known clinicopathological risk factors such as the EORTC score. Univariate Cox analysis identified age and tumor stage, along with the four risk estimators, as independent prognostic indicators of PFS in NMIBC (Table 5). Because tumor stage is an index used to calculate the EORTC, we excluded the EORTC from the following multivariate Cox analyses. Multivariate analysis based on tumor age, stage and the MSP888 classifier, or the UROMOL classification, or the Lund subtype revealed that the MSP888 classifier and the Lund subtype retained statistical significance with respect to PFS of NMIBC patients (HR, 3.025 and 3.605; 95% CI, 1.277-7.166 and 1.200-10.826; and p = 0.012 and 0.022, respectively) ( Table 5). Kaplan-Meier analysis revealed that NMIBC patients in Cluster 3 of MSP888, in class 2 of the UROMOL classification or with a genomically unstable Lund subtype experienced more progression than patients classified as other subtypes (log-rank test, p < 0.05) (Figure 3), confirming the high prognostic relevance of the three clusters classified by the MSP888.    Finally, we evaluated the prognostic significance of the four risk predictors for NMIBC progression in patients not treated with BCG. Univariate and multivariate Cox regression analyses revealed that MSP888 classifier and Lund subtype were independent predictors of PFS. The risk of progression for patients in Cluster 3 of the MSP888 was 3.753 times higher than that for patients in Clusters 1 and 2 (p = 0.005). In addition, the progression risk for patients with a genomically unstable Lund subtype was 5.962 times higher than that of patients with other subtypes (p = 0.005) ( Table 6). The PFS of NMIBC patients not treated with BCG was estimated from Kaplan-Meier survival plots. Patients in Cluster 3 of the MSP888, in class 2 of UROMOL or with a genomically unstable Lund subtype experienced more progression than other subgroups (log-rank test, p < 0.05) (Figure 4), suggesting the crucial prognostic relevance of the three clusters classified by the MSP888, which is comparable with that of previously constructed global subtypes.

Discussion
Most traditional classification systems for BCa are based on pathological parameters [15,16]. However, there are large differences between individuals with respect to prognosis of BCa, even in those with tumors of similar pathological stage and grade [17,18]. Accordingly, an omnibus risk score table (i.e., EORTC), which is based on clinicopathologic characteristics, was developed in 2006 to assess the probability of recurrence and progression in patients with NMIBC [5]. However, the lack of patients treated with BCG led to development of the CUETO risk tables to predict the risk of recurrence and progression specifically in patients treated with BCG [6]. To estimate 1−year and 5−year recurrence or

Discussion
Most traditional classification systems for BCa are based on pathological parameters [15,16]. However, there are large differences between individuals with respect to prognosis of BCa, even in those with tumors of similar pathological stage and grade [17,18]. Accordingly, an omnibus risk score table (i.e., EORTC), which is based on clinicopathologic characteristics, was developed in 2006 to assess the probability of recurrence and progression in patients with NMIBC [5]. However, the lack of patients treated with BCG led to development of the CUETO risk tables to predict the risk of recurrence and progression specifically in patients treated with BCG [6]. To estimate 1-year and 5-year recurrence or progression rates, the EORTC and CUETO risk scores classify patients into four groups (low, intermediate low, intermediate high and high). Several external validation studies have been performed to evaluate the applicability of EORTC and CUETO in different countries, although the results are controversial [7]. For example, some patients were pathologically classified as low risk, but the actual tumor showed highly invasive biological characteristics and early metastasis. Therefore, physicians can find it difficult to provide individualized treatment and management strategies for BCa patients. Nevertheless, to date, these two scores are the most common and best-validated models for predicting BCa prognosis. In late 2021, EAU presented a restructured risk group that successfully stratified progression risks in NMIBC [8]. Despite the limitation of overestimating progression in patients with BCG treatment, this new risk stratification is worth looking forward to.
Analyses of genomic alterations in NMIBC revealed complex genomic patterns underlying bladder carcinogenesis; however, the pathological parameters of the tumor do not fully reflect the "intrinsic characteristics" [12,15,17,18]. Earlier studies on BCa performed transcriptomic analyses to develop molecular classification systems for NMIBC [9,12,13,19,20]. The UROMOL project is the most well-known study of NMIBC; the aim was to predict the disease course of BCa using risk scores that combine molecular and clinical risk factors. The UROMOL classification was developed in 2016 as it was based on data from 460 NMIBC patients. The UROMOL 2016 classification system has three gene expression-based classes: classes 1-3. Each class has a different clinical outcome and molecular characteristics [9]. A more recent study reported an integrative multi-omics analysis of NMIBC from a total of 834 patients included in the UROMOL project, and identified four classes (1, 2a, 2b and 3) that reflect tumor biology and disease aggressiveness [12]. Although previous studies identified the molecular characteristics of NMIBC, they provided limited information concerning the association between molecular taxonomy and therapeutic relevance (particularly responses to BCG treatment). Previously, we developed and validated a MSP888 predictor for classifying prognostic subtypes of NMIBC through gene expression profiling of multiple cohorts comprising 948 NMIBC patients [13]. Three subtypes of NMIBC exhibited distinct prognostic features: (1) DP.BCG+ was associated with progression and a positive responsive to BCG; (2) REC.BCG+ was associated with worse RFS and a better response to BCG; and (3) EP exhibited equivocal prognostic behavior.
Intravesical BCG therapy is the most recommended option for intermediate-and highrisk NMIBC patients after transurethral resection of bladder tumor (TURBT) to reduce the risk of NMIBC recurrence and progression [1,3,21]. With BCG treatment, approximately 35% of patients could experience long-term, sustained remissions, although half of patients will suffer recurrence within two years, and about one-third will progress to muscleinvasive stages ultimately requiring radical cystectomy to remove the bladder [3,21]. Risk stratification usually relies on the clinicopathologic features of the disease to provide an ideal therapy for each individual patient [1,3,21]. However, this stratification often underestimates those patients assessed as low-risk NMIBC who were not recommend the BCG treatment. The right risk estimation leads to the right disease management, which directly correlates with the patient's prognosis and their overall survival. In order to more effectively tailor an appropriate treatment regimen for these low-or intermediate-low-risk patients, a more reliable further stratification is needed [21].
In the current study, we compared the predictive value of the MSP888 with that of the EORTC, CUETO and 2021 EAU risk scores for NMIBC recurrence or progression in a CBNUH cohort (n = 178). Furthermore, we compared the predictive value of the MSP888, UROMOL classification, Lund subtype and EORTC risk score for NMIBC progression in the UROMOL cohort (n = 460). The results from the CBNUH cohort demonstrated that the MSP888 classifier is an independent risk factor for NMIBC prognosis, and that its per-formance is superior to that of the EORTC, CUETO and 2021 EAU risk scores, especially for NMIBC patients not treated with BCG (Tables 1-3, Figure 2). In the UROMOL cohort, the ability of the MSP888 classifier to predict NMIBC progression was comparable with that of the EORTC score, as well as with the UROMOL classification and Lund subtype (Table 5, Figure 3). In particular, the MSP888 classifier and Lund subtype were independent risk factors for NMIBC progression in patients not treated with BCG (Table 6, Figure 4). Kaplan-Meier survival plots revealed that all four risk estimators (MSP888, UROMOL classification, Lund subtype and the EORTC score) could discriminate NMIBC patients with progression from those without progression (log-rank test, p < 0.05 respectively); however, the EORTC result was equivocal. The EORTC scores divided patients into two groups, which is in contrast to the common division into four groups; in addition, the number of tumors (an index used to calculate the EORTC score) was not provided for the UROMOL cohort. Accordingly, these results suggest that our molecular signature-based risk predictor MSP888 could provide a more accurate risk assessment than clinicopathological-based risk scores for those patients not treated with BCG, most of whom were evaluated as lowor intermediate-low-risk NMIBC.
The current study is a second phase validation of the MSP888 predictor for NMIBC patients. However, the small sample size of RNA-seq is the limitation of the present study. Furthermore, a series of tests will be continuously arranged to verify the predictive power of MSP888. Finally, we aimed to manufacture a prognostic molecular subtypebased panel based on counsel to NMIBC patients about their risk of recurrence or progression, followed by the suggestion of individual treatment strategies. In the next generation of risk classification, we are likely to see the inclusion of molecular subtyping with specific treatment considerations.

Study Design
Overall study design is shown in Figure 5. To determine the prognostic utility of the signature-based classifier MSP888, RNA-seq data from 49 NMIBC patients registered at CBNUH were used to establish the MSP888 classifier. The performance of the MSP888 classifier for predicting NMIBC recurrence and progression in 178 NMIBC patients from CBNUH (including 49 patients with newly analyzed RNA-seq data and a published cohort of 129 that was used to develop the subtype classifier and to validate its prognostic value [13]) was compared to that of EORTC, CUETO and 2021 EAU risk groups. In addition, the prognostic ability of the MSP888 classifier, EORTC, the UROMOL classification and Lund subtype were evaluated in 460 NMIBC patients from the UROMOL cohort (a European multicenter prospective study) [9]. In particular, the risk for patients with and without BCG treatment in both the CBNUH and UROMOL cohorts was estimated.

Patients and Tissue Samples
The study methods conformed with the standards set out by the Declaration of Helsinki. Table 7 shows the baseline characteristics of the NMIBC subjects used for RNA-seq analysis. Forty-nine NMIBC patients with distinct clinical outcomes were selected. Tumor tissues were collected from surgically-resected NMIBC from September 2001 to March 2018. All tumors were macro-dissected, typically within 15 min of surgical resection. Each NMIBC specimen was confirmed by pathological analysis of a part of a fresh-frozen specimen obtained from TURBT, and was histologically verified as transitional cell carcinoma. The biospecimens were obtained from the CBNUH (Cheongju, Republic of Korea), which is involved in the National Biobank of Korea. The study was approved by the Institutional Review Board at CBNUH (GR2010-12-010 and GR2020-07-018), and the experiments were undertaken with the informed written consent of all participants. To reduce the chances of confounding factors affecting the analyses, patients diagnosed with concomitant CIS, or CIS lesions alone, were excluded. Tumors were staged (2017 TNM Classification) and graded (1997 WHO Classification) according to standard criteria [3]. The EORTC and CUETO risk scores for both recurrence and progression were calculated [5,6], and patients were stratified into four risk groups according to the calculated scores. Each patient was followed up and managed according to standard guideline recommendations [3]. Surveillance was performed by cystoscopy and upper urinary tract imaging in accordance with European Association of Urology guidelines [3]. In the current investigation, disease recurrence was defined as relapse of primary NMIBC of the same pathologic stage, and progression of NMIBC was defined as TNM stage progression after disease recurrence [22]. The mean follow-up period for NMIBC patients was 64.43 months (range, 4.90-227.33).
with specific treatment considerations.

Study Design
Overall study design is shown in Figure 5. To determine the prognostic utility of the signature−based classifier MSP888, RNA−seq data from 49 NMIBC patients registered at CBNUH were used to establish the MSP888 classifier. The performance of the MSP888 classifier for predicting NMIBC recurrence and progression in 178 NMIBC patients from CBNUH (including 49 patients with newly analyzed RNA−seq data and a published cohort of 129 that was used to develop the subtype classifier and to validate its prognostic value [13]) was compared to that of EORTC, CUETO and 2021 EAU risk groups. In addition, the prognostic ability of the MSP888 classifier, EORTC, the UROMOL classification and Lund subtype were evaluated in 460 NMIBC patients from the UROMOL cohort (a European multicenter prospective study) [9]. In particular, the risk for patients with and without BCG treatment in both the CBNUH and UROMOL cohorts was estimated.

Patients and Tissue Samples
The study methods conformed with the standards set out by the Declaration of Helsinki. Table 7 shows the baseline characteristics of the NMIBC subjects used for RNA−seq analysis. Forty−nine NMIBC patients with distinct clinical outcomes were selected. Tumor tissues were collected from surgically−resected NMIBC from September 2001 to March 2018. All tumors were macro−dissected, typically within 15 min of surgical resection. Each NMIBC specimen was confirmed by pathological analysis of a part of a fresh−frozen specimen obtained from TURBT, and was histologically verified as transitional cell carcinoma. The biospecimens were obtained from the CBNUH (Cheongju, Republic of Korea), which

RNA Extraction
Total RNA was extracted from tissues using TRIzol reagent (Invitrogen, Carlsbad, CA, USA), as described previously [23], and stored at −80 until use.

RNA-Sequence Analysis
The value of RNA Integrity Number (RIN) and the DV200 metric were measured using an RNA 6000 Nano Kit and an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) to confirm the quality and integrity of the RNA. RNA samples with a RIN higher than seven were designated as "good total RNA quality" and selected for downstream application. Total RNA (500 ng) was processed to prepare a whole transcriptome sequencing library. Whole transcriptome RNA was enriched by depleting ribosomal RNA (rRNA) prior to generating the whole transcriptome sequencing library using the MGIEasy RNA Directional Library Prep Kit (MGI Tech Co., Ltd., Shenzhen, China). After rRNA was depleted, the remaining RNA was fragmented into small pieces using divalent cations under elevated temperature. The cleaved RNA fragments were copied to generate first strand cDNA using reverse transcriptase and random primers. Strand-specificity was achieved using an RT directional buffer. This was followed by second strand cDNA synthesis. The cDNA fragments had a single "A" base added prior to ligation of the adapter. The products were then purified and enriched by PCR to create the final cDNA library. The doublestranded library was quantified using a QauntiFluor ONE dsDNA System (Promega, Madison, WI, USA) and add TE buffer to equal 330 ng in a total volume of 60 . The library was cyclized at 37 • C for 60 min and then digested at 37 • C for 30 min, followed by cleanup of the circularization product. To generate a DNA nanoball (DNB), the library was incubated at 30 • C for 25 min with DNB enzyme. Finally, the library was quantified using a Qaunti-Fluor ssDNA System (Promega, Madison, WI, USA). Sequencing of the prepared DNB was conducted using the MGIseq system (MGI Tech Co., Ltd., Shenzhen, China) with 150 bp paired-end reads. Reference genome sequence data from Homo sapiens were obtained from the NCBI Genome database (assembly ID: GRCh38). Reference genome indexing and read mapping of tissue samples were performed using STAR software (ver. 2.5.4b) [24].

Transcriptomic Summarizing
All of the gene expression data were log2-transformed and quantile-normalized in the current study. The read counts per million fragments mapped value of each sample were measured in the 49 NMIBC patients of the CBNUH-RNAseq cohort. Mean sequencing coverage of RNA-seq data was 47.27× ± standard deviation 16.39× (ranges from 32.70× to 95.10×) and the average mapping rate to reference genome is 93.07% ± standard deviation 1.46% (ranges from 88.86 to 95.20%).

Establishment of MSP888 Clusters in the CBNUH RNA-seq Cohort
In a previous study [13], we generated the MSP888 prediction model by adopting a DBN [14] algorithm. MSP888 is a model incorporating 888 differentially expressed genes among three subtypes according to two sample t-tests. Then, 5 hidden layers were set up, in which 600, 300, 100, 300 and 600 nodes were allocated respectively to construct a fully connected neural network. A Tanh function was used as an activation function for the neural network, and a training procedure was repeated during 1000 epochs. To prevent overfitting or non-convergence of the classifier, the prediction model was pre-trained by an AutoEncoder [14] algorithm. Development of the prediction model was undertaken using the H2O (https://www.h2o.ai (accessed on 3 December 2021)) deep learning platform (ver. 3.32.0.2). Using this training model, evaluation of predicted outcomes in the CBNUH RNAseq cohort was carried out based on gene expression signatures.

Public Datasets of NMIBC Patients
An RNA-seq dataset from the European UROMOL consortium [9] was downloaded from the ArrayExpress database under accession number E-MTAB-4321. In this dataset, 460 NMIBC samples were used as a validation cohort (the UROMOL cohort). All transcriptomic data used in this study include patient survival and follow-up time to estimate the prognostic relevance of the signatures.

Statistical Analyses
The relationship between each risk predictor and NMBC prognosis was analyzed using Pearson's Chi-squared test. The interval LR for each stratum was calculated as the likelihood of that test result in patients with a positive test divided by the likelihood of that result in patients with a negative test; this was done to show how likely it is that a patient will experience recurrence or progression. The significance of various clinicopathological variables and the risk estimators was evaluated using univariate and multivariate Cox proportional hazard regression models. Relative risk was determined by calculating hazard ratios (HRs) and 95% confidence intervals (CIs). Kaplan-Meier survival curves were plotted to examine the prognostic value of each risk estimator and compared using the logrank test. Statistical analyses were performed using IBM SPSS Statistics ver. 24.0 (IBM Co., Armonk, NY, USA) and MedCalc ver. 18.2.1 (MedCalc Software, Mariakerke, Belgium). p-Values < 0.05 were considered significant.

Conclusions
In conclusion, the current study validated the prognostic power of the MSP888 classifier and compared it with that of the EORTC, CUETO and 2021 EAU risk scores. The MSP888 classifier showed superior value for predicting NMIBC prognosis, suggesting that molecular classification-based systems are more accurate than the risk scores based on clinicopathological characteristics. This explains the fact that, compared to the traditional classification system, molecular subtyping reflects the intrinsic characteristics of the tumors, and predicts the prognosis and treatment response of NMIBCs.

Institutional Review Board Statement:
The study was conducted in accordance with the Declaration of Helsinki and approved by the Institutional Review Board at Chungbuk National University (GR2010-12-010 and GR2020-07-018). All samples derived from the National Biobank of Korea were obtained with informed consent under Institutional Review Board-approved protocols.
Informed Consent Statement: All samples derived from the National Biobank of Korea were obtained with informed consent under Institutional Review Board-approved protocols.

Data Availability Statement:
The datasets used and/or analyzed in the current study are available from the corresponding authors upon reasonable request.