Proteogenomic characterization of the non-muscle-invasive bladder cancer response to BCG reveals potential therapeutic strategies

doi:10.21203/rs.3.rs-4008035/v1

Download PDF

Research Article

Proteogenomic characterization of the non-muscle-invasive bladder cancer response to BCG reveals potential therapeutic strategies

https://doi.org/10.21203/rs.3.rs-4008035/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background

Intravesical bacillus Calmette-Guérin (BCG) is the standard therapy for adjuvant treatment in patients with intermediate- and high-risk superficial bladder cancer. However, the molecular properties associated with BCG therapy have not been fully characterized.

Methods

We reported a comprehensive proteogenomic analysis, including whole-genome sequencing, proteomics, and phosphoproteomics profiling, of 160 non-invasive-muscle bladder cancer (NMIBC) patients treated with BCG.

Results

Proteogenomic integration analysis indicted that tumor mutational burden (TMB), associated with STAT1 activity, was relevant to drug sensitivity. Additionally, our analysis of copy number alterations (CNAs) showed that TLR3 deletion was negatively correlated with response to BCG therapy. TLR3 was validated to regulate the cytokine secretion, and enhance sensitivity to BCG in BC cell lines and organoids. High TMB levels were also associated with improved BCG efficacy across different TLR3 expression subgroups, which holds significant implications. Through proteomic analysis, we identified three subtypes in patients with BCG, reflecting distinct clinical prognosis and biological characteristics. Furthermore, we developed prognostic models with high accuracy to predict the therapeutic response and PFS of NMIBC.

Conclusions

This study provides a rich resource for investigating the mechanisms and indicators of BCG therapy in NMIBC, which can be basis for further improvement of therapeutic response.

Bladder cancer (BC) is the tenth most commonly diagnosed cancer worldwide [1]. There are approximately 570,000 new cases and 210,000 related deaths reported in 2020 [2]. The most common presenting symptom observed in patients with bladder cancer is microscopic or gross hematuria, and risk factors for developing bladder cancer consists of male sex, smoking, exposure to certain drug, and diabetes [3, 4]. The most prevalent pathological type of BC is urothelial carcinoma (UC), with approximately 75% of patients diagnosed with non-muscle-invasive bladder cancer (NMIBC). NMIBC includes Ta, T1, and carcinoma in situ (CIS). While the 5-year overall survival (OS) rate for NMIBC is 90%, the recurrence rate ranges from 60–70% [5].

Treatment strategies for NMIBC focus on reducing recurrence and preventing disease progression. According to the risk stratification by the American Urological Association (AUA) and the Society of Urologic Oncology (SUO), NMIBC can be categorized into three groups: Low- risk, Intermediate-risk, and High-risk patients. Although the survival rate is favorable, patients with Low- and Intermediate-risk NMIBC have a 5-year recurrence-free survival rate of 43% and 33%, respectively, and up to 21% with High-risk disease progress to MIBC [6, 7]. Bacillus Calmette-Guerin (BCG), a treatment option used since the 1970s, could reduce the risk of recurrence and progression in BC [8]. For Intermediate-risk and high-risk NMIBC patients, Intravesical BCG is the first choice. A meta-analysis of NMIBC studies demonstrated that BCG improved the recurrence-free rate from 26–47% (odds ratio, 0.41[95% Cl, 0.30–0.56]) and the progression-free (PFS) rate from 20–15% (OR, 0.74 [95% Cl, 0.45–1.22]) compared to various intravesical chemotherapy agents [9].

Although BCG immunotherapy has been used to treat NMIBC for more than 40 years, the mechanisms by which BCG immunotherapy mediates tumor immunity remains unclear [10, 11]. Recent studies have suggested that BCG-induced specific tumor immunity seems likely because animals rechallenged with tumor cells following tumor challenge and BCG therapy mount an adaptive response and effectively eliminate these tumor cells [12, 13]. However, there is currently no clear data available in humans. Deeper understanding of the mechanisms behind BCG’s antitumor activity is essential to improve treatment for NMIBC.

Approximately 20–50% of patients, depending on follow-up time and risk profile, experience a recurrence after receiving adequate BCG therapy [14, 15]. Biomarkers that can predict the response to BCG in NMIBC might benefit both patients and clinicians in making treatment-related decisions. A model that can predict the risk of recurrence and progression after 12 doses of intravesical BCG has been published by the Club Urologico Espanol de Tratamiento Oncologico (CUETO). Using this model, which include tumor size, tumor focality, grade and tumor stage, the calculated risk of recurrence is lower than that obtained from the EORTC risk tables [16]. However, the EORTC risk tables are still not considered valid for use in patients alone. It is important to note that these models have their own limitations. Better markers of immune response and molecular profiling of BCG therapy are urgently needed to prevent ineffective therapy and enable a personalized therapy.

Several molecular classifications of NMIBC have been proposed to establish a connection between molecular patterns and clinical features. Lindskrog and his colleagues identified four classes (1, 2a, 2b, and 3) reflecting tumor biology and disease aggressiveness in a study involving 834 NMIBC patients. They observed that class 2b tumors, with the highest PD-L1 expression, may respond poorly to BCG treatment [17]. In another study, Robertson and his colleagues performed transcriptome profiling on 73 high-risk NMIBC patients treated with BCG and identified five consensus subtypes. They further found that T1-Myc and T1-Early subtypes had the highest recurrence rates, with over half of the tumors recurring after BCG treatment. [18]. However, these classifications are mainly based on transcriptional data, and there is limited research on proteome-based classifications for NMIBC treated with BCG. With proteins being directly linked to phenotypes, protein-based molecular subtyping shows promise in providing critical information on translating genome signals to cell function. In contrast to the extensive research on transcriptomics, deep proteomic analysis of larger cohorts treated with BCG are need to delineate response mechanisms. Additionally, systematic data regarding the outcomes of NMIBC patients undergoing BCG therapy has not been reported.

In this study, we reported an integrative multi-omics analysis of NMIBC tumors from 86 patients treated with BCG. Proteogenomic integration analysis indicated that high levels of TMB and TLR3 were associated with a good prognosis, through a common mechanism of modulating antigen processing and presentation. In addition, proteomic analysis identified three subtypes reflecting distinct prognostic and biological characteristics. These subtypes exhibited different protein signatures and kinase activity which were found to be associated with differential responses to BCG therapy. Furthermore, we developed prognostic models with high accuracy to predict the therapeutic response to BCG and the PFS of NMIBC. Collectively, this study portrayed a NMIBC cohort of proteogenomic landscape and provided valuable insights into BCG therapy, which could help in better management of NMIBC.

Proteogenomic characterization of NMIBC

To obtain a comprehensive molecular understanding of Chinese NMIBC patients treated with BCG, frozen tumors from 160 NMIBC patients (discovery cohort [n = 86], validation cohort [n = 74]) were selected for proteogenomic analysis based on stringent criteria ((Fig. 1A; Table S1). These tumors were collected from newly diagnosed NMIBC patients who underwent surgical resection without any prior treatment such as chemotherapy or radiotherapy. Following transurethral resection of bladder tumor (TURBT), all patients received BCG instillations at standard dosing as described in previous studies [19]. Clinical follow-up was conducted according to the guidelines of the European Associated of Urology (EAU), and detailed information on resources and patient inclusion can be found in the Methods section. All samples were histologically scored by two expert genitourinary pathologists according to the EAU guidelines on NMIBC and the samples were grouped into response (R) and non-response (NR) categories. Among the 86 NMIBC patients in discovery cohort, the recurrence rate was 45.3% (39/86), with progression observed in 4.7% (4/86). A schematic of the experimental design is depicted in Fig. 1A-B. The baseline characteristics of the groups are shown in Table 1, and no differences in age, gender, or tumor stage at initial TURBT were identified between the BCG-R and BCG-NR groups.

Proteomics measurement of all patient samples resulted in a total of 14,956 at a 1% false discovery rate (FDR) at the protein and peptide levels (Figure S1A). 8,389, 7,982, and 6,298 protein groups were identified in BCG-R, BCG-NR, and normal samples, respectively (Figure S1B). The tumor proteome and the normal proteome exhibited a unimodal distribution (Figure S1D). Furthermore, a total of 92,524 phosphosites corresponding to 8,972 phosphoproteins were identified (Figure S1C). The number of genes identified as corresponding to the proteome and phosphoproteome was 8,881 (Figure S1E). The average Spearman’s correlation coefficient, calculated for all quality control runs of HEK293T cell samples, was 0.9, indicating high quality of the mass spectrometry (MS) data (Figure S1F). Additionally, the average correlation coefficient in the phospho-proteome was 0.92 (Figure S1G), indicating the stability of our MS platforms.

Using whole-exome [20] sequencing data, we identified 17,925 genetic variation events including 11,936 single-nucleotide variants (SNVs) and 9,765 small insertion-deletions (indels) through the comparison between tumor and blood (Fig. 1C). Most of the previously reported genomic findings for NMIBC were confirmed. The two most frequently mutated chromatin-modifying genes are KMT2D and KDM6A. Some well-known bladder cancer-related genes, such as TP53, FGFR3, PIK3CA and ERCC2, were also detected in our cohort (Fig. 1C). Additionally, the frequency of these hotspot mutations was similar to what has been previously reported [21], including KDM6A (38%), KMT2D (34%), ERCC2 (22%), and ATM (13%). Furthermore, to identify the genomic alterations associated with recurrence after BCG, we compared the mutational frequencies of hotspot mutations between BCG-R and BCG-NR, the result showed no significance (Fig. 1D).

The presence or absence of a mutational signature can aid in the clinical interpretation of potentially actionable mutations. Because mutational signatures describe the state of the DNA repair machinery of cancer cells, they can potentially serve as drug sensitivity marker [22]. We further identified four mutational signatures by using Sigminer (Fig. 1E). These signatures, namely Signature 1–4 corresponded to the known COSMIC (Catalog of Somatic Mutations in Cancer) Signatures: SBS2, SBS13, SBS5, SBS12. We compared the tumor mutation burden (TMB) associated with these four mutational signatures and found that SBS2, which represents the activity of APOBEC, had the highest TMB level (Fig. 1F). In addition, survival analysis showed that SBS2 was associated with the best prognosis (Fig. 1G). These findings suggested that NMIBC patients with high TMB have a better prognosis.

High mutational burden (TMB) is relevant to drug sensitivity

BC has a high rate of nucleotide variants (Figure S2A), indicating that it might benefit from immunotherapy. TMB is emerging as an important factor in cancer immunotherapy as it predicts responses in various types of cancer [23]. In this study, we investigated the association between TMB and clinical benefit from BCG immunotherapy in NMIBCs. The results showed that the TMB value was higher in BCG-R compared to BCG-NR (Fig. 2A, Kruskal–Wallis test, p < 0.05). While there were no significant differences in WES coverage depth and tumor purity between patients with BCG-R and BCG-NR (Figure S2B and S2C). We performed a stratified analysis by selecting the highest quartile of mutation (top 25%) as the TMB-high group (defined as equal to or above the median TMB of the cohort: ≥10 mutations/Mb). Using this approach, we observed that BCG-R was higher in patients with high TMB, which was significantly associated with response to BCG immunotherapy, compared to low TMB (Fig. 2B, Fisher’s exact test). We further found that NMIBC patients with high TMB had an improved PFS compared to those with low TMB (Fig. 2C, log-rank test, p < 0.05). Additionally, similar stratified analyses based on TMB were performed in the Fudan cohort (Figure S2D, log-rank test, p < 0.05) [24] and the MSKCC cohort (Figure S2E, log-rank test, p < 0.05) [23], both of which showed that patients with high TMB had a better prognosis.

Moreover, the prediction accuracy of BCG response was calculated based on clinical features, such as tumor stage, tumor size, and tumor number, either individually or in combination with TMB. The area under the curve (AUC) of the combination of clinical features and TMB was found to be higher than that of clinical features alone (Fig. 2D). These results indicated that BCG therapy may be more beneficial for patients with higher TMB, and TMB could serve as a potential predictor of BCG responses.

To investigate the potential association between TMB and improved survival after BCG therapy, we tried to establish a connection between TMB and its clinical features. Firstly, we found that TMB was positively correlated with DNA replication, mismatch repair, and cell cycle pathways (Fig. 2E). Additionally, when we focused on immune-related pathways, we observed a significant positive correlation between TMB and the JAK STAT signaling pathway. Furthermore, TMB was observed to be positively associated with predicted STAT1 activity, which was inferred based on the expression of its target genes (Fig. 2F, Spearman correlation r = 0.38, p = 0.0011; Methods). It has been reported that STATs play a complicated and fundamental role in the immune system, and have been closely linked to immunotherapy response [25, 26]. Furthermore, the STAT1 activity was higher in TMB high group (Fig. 2G, Wilcoxon rank-sum test, p < 0.05), and this increased activity was associated with poor prognoses (Fig. 2H, log rank test, p < 0.05).

To further explore the relationship between STAT1 activity and BCG response, we performed analysis on the targets of STAT1. We observed that the expression of many proteins, which are known to be STAT1 targets increased as STAT1 activity increased (Fig. 2I). Proteasome 20S subunit beta 9 (PSMB9), one of these STAT1 targets, is an interferon-gamma (IFN-γ) inducible factor located within the major histocompatibility complex (MHC) class II region on chromosome 6 [27]. PSMB9 is directly regulated by the anti-oncogenic IRF1, a well-known gene with tumor suppressor activity [28–30]. Notably, ssGSEA analysis indicated PSMB9 was significantly positive with immune-related pathways, such as antigen processing and presentation, proteasome, and chemokine signaling pathway (Fig. 2J). Further analysis showed that the proteins of the proteasome were upregulated in conjunction with TMB (Fig. 2K and 2L). These results indicated that high STAT1 activity might upregulate antigen processing and presentation though its target, PSMB9, which improved the therapeutic efficacy of BCG.

Taken together, the TMB value in BCG-R was higher in BCG-NR, and high TMB was associated with a better outcome in NMIBC patients who received BCG therapy. Further analysis showed that high TMB was correlated with increased STAT1 activity, which might enhance antigen processing and presentation by PSMB9.

The expression level of TLR3 in tumors is associated with prognosis and BCG treatment response

Somatic copy-number alterations (SCNAs) based on WES data showed the most frequent gains in chromosomes 1q, 8q, 20p, and 20q as well as losses in chromosomes 9p, 9q, and 17p at the arm level of SCNAs (Figures S3A; Table S3). Focal level SCNAs included amplifications of 1q23.3, 1q21.2, 3q27.1, and 11q13.3 and deletions of 4q35.1, 5q13.1, and 9p21.3, among others (Fig. 3A and 3B; Table S3). These SCNAs events are consistent with previous reports in the NMIBC cohort [17, 31]. We found that amplifications of 1q23.3, 1q21.2 and deletion of 4q35.1 were associated with poor patient prognoses (Fig. 3C). We further examined the regulatory effects of these focal level SCNAs on protein expression. SCNAs can affect protein abundance through either “cis” or “trans” modes. We observed significant positive cis-effects on protein levels for OTUD7B on peak 1q21.2, F11R on peak 1q23.3, and TLR3 on peak 4q35.1, among others (Fig. 3D). However, only TLR3 was associated with prognoses among these nine cis-effects (Fig. 3E). Patients with higher expression of TLR3 protein in tumors appeared to have better prognostic outcomes (Fig. 3F, log rank test, p < 0.05). This association was confirmed in the Chicago T1BC cohort [18], which was a BCG treatment cohort as well (Figure S3B). TLR3 is predominantly expressed in innate immune cells such as macrophages, dendritic cells (DCs) and natural killer (NK) cells and is involved in the activation of innate immunity, which results in the production of proinflammatory cytokines, chemokines, and adhesion molecules [32]. To further confirm the expression of TLR3 in NMIBC tumors, immunohistochemistry (IHC) staining was performed on tumor tissues. Positive staining for TLR3 was observed in NMIBC tumors (Figure S3C), consistent with previous findings [33]. In addition, we found that TLR3 was overrepresented in BCG-R tumors compared to BCG-NR tumors (Figure S3C), which indicates that detecting TLR3 expression may help identify patients who respond to BCG.

To better understand the effect of TLR3 on prognosis and BCG treatment response, we explored the correlation between the protein abundance of TLR3 and pathways. The abundance of TLR3 protein in tumors was positively correlated with the antigen processing and presentation gene set in the KEGG (Kyoto Encyclopedia of Genes and Genomes) database (Fig. 3G; Table S3). As the TLR3 protein abundance increased, the expression of human leukocyte antigen (HLA) proteins from the MHC class I or II also increased (Figure S3D). Furthermore, TLR3 is normally located in endosomes and binds to ligands, leading to the activation of transcription factors [34, 35]. Therefore, we further focused on the correlation between the protein abundance of TLR3 and transcription factors. NF-κB1, an important transcription factor associated with TLR3 signaling [34], showed the highest positive correlation with the abundance of TLR3 protein in tumors (Fig. 3H; Spearman’s r = 0.32, p = 7.3 × 10E-4). To confirm the correlation between TLR3 and NF-κB1, we enrolled the TCGA BLCA cohort data and examined the association between TLR3 and NF-κB1. Consistently, a significantly positive correlation between NF-κB1 and TLR3 was also identified in the TCGA BLCA cohort (Figure S3E; Spearman’s r = 0.57, p = 0), indicating an interaction between TLR3 and NF-κB1 in urothelial bladder cancer. Moreover, patients with higher expression of NF-κB1 protein in tumors appeared to have better prognostic outcomes (Figure S3F; log rank test, p < 0.05). NF-κB is responsible for physiological functions like inflammation, proliferation, survival, and stress responses [36]. We also found that the abundance of NF-κB1 protein in tumors was positively correlated with the pattern recognition receptor signaling (RIG, NOD, and TOLL like) as well as antigen processing and presentation gene sets (Fig. 3I and 3J). Conversely, the gene sets related to the multi-gene proliferation score (MGPS; Methods) showed a negative correlation (Figure S3G).

The biological functioning of NF-κB requires dimerization, activation and translocation to the nucleus for gene expression [37]. The phosphorylation of IκBs by IKKs is a crucial step in the activation of NF-κB [38]. Notably, positive associations between IKBKB and IKBA at Ser 32 and Ser 36 were observed in our phosphoproteome data (Fig. 3K). The expression of many proteins involved in the antigen processing and presentation, which are targets of NF-κB, showed favorable prognostic signatures and upregulation with increasing TLR3 protein abundance (Fig. 3J). In addition, elevated antigen processing and presentation pathway signaling was associated with increased survival (Figure S3H; log rank test, p < 0.05). Since we also observed an association between high TMB and BCG response (in Fig. 2), we investigated whether there is a synergistic effect between TMB and TLR3. Interestingly, we found that patients with high levels of both TMB and TLR3 had better survival than those without high levels in tumors (Fig. 3L; log rank test, p < 0.05). These results suggested that both high levels of TMB and TLR3 contribute to a good prognosis, possibly through a common mechanism of modulating antigen processing and presentation (Fig. 3M).

Proteomic subtyping of the NMIBC BCG cohort and their association with the therapeutic response

Consensus clustering identified three proteomic subtypes based on the 2,348 most variable proteins (proteins with the top 25% standard deviations) that were present in more than 50% of 73 pre-BCG tumor samples (Fig. 4A and S4A). They were designated as GP1 (n = 29), GP2 (n = 35), and GP3 (n = 9), and they exhibited distinct molecular and clinical features (Fig. 4A; Table S4). Patients in GP1 and GP2 had better clinical outcomes, whereas patients in GP3 had the worst prognosis and the highest percentage of BCG non-responders (Fig. 4B-4D; log rank test, p < 0.05). Combined with the clinical data, tumors with papillary were mostly enriched in the proteomic subtype GP2 (Fig. 4A; Fisher’s exact test, p < 0.05), but no association was observed between proteomic subtypes and tumor grade, tumor stage, gender, or concomitant CIS (p > 0.05). Notably, the proteomic subtype GP3 exhibited the lowest levels of TMB and TLR3 (Fig. 4A and 4E), which was consistent with the observation that high levels of TMB and TLR3 were associated with a good prognosis. To explore the overlap of the proteomic subtypes with previously published subtypes, we compared them with the Chicago T1BC subtypes [39], UROMOL NMIBC subtypes [40], UNC MIBC subtypes [41], and Consensus MIBC subtypes [42] (Fig. 4A and S4B; Table S4). In the Chicago T1BC subtypes, subtype GP3 overlapped with the easily recurring Early subtype, while subtype GP1 overlapped with the LumGU subtype, which has a better prognosis (Fig. 4A and S4B, Fisher's exact test, p < 0.05). In the UNC MIBC subtype and Consensus MIBC subtypes, subtype GP1/2 were predominantly luminal or luminal papillary subtypes, which are known to have a better prognosis. These results indicated that our proteomic subtypes exhibit the similar survival patterns when compared to other subtypes.

Subtype-specific pathway enrichment analysis indicated different features among the three proteomic subgroups. Subtype GP1 was characterized by the highest level of immune-related pathways, such as chemokine signaling, antigen processing and presentation, and complement and coagulation cascades (NCF1, FGG, PLCG2, etc.) (Fig. 4F and 4G). On the other hand, subtype GP2 was more closely associated with transcription, including transcriptional mis-regulation in cancer, transcription initiation from RNA polymerase II promoter (MECOM, TAF5, SIN3A, etc.), and inositol phosphate metabolism. Subtype GP3 was characterized by translation and ribosome biogenesis pathways (RPS16, RPL27, LTV1, etc.), such as ribosomal small subunit biogenesis, translation, and mRNA processing pathway (Fig. 4F and 4G). In addition, the expression levels of several HLA proteins from MHC class I and II were found to be the lowest in in the subtype GP3 (Figure S4C; Table S4). These results suggested a possible mechanism by which cancer cell evade the immune system and influence the outcome of BCG immunotherapy.

Comparison of the phosphoproteome among the three proteomic subtypes revealed subtype-specific activated kinases (Methods), including PRKACA and FGF for GP1, ALB2 and HIPK2 for GP2, and PASK and CLK4 for GP3 (Figure S4D). Further investigation into the differentially altered phosphosites showed that the phosphorylation of the PASK and CLK4 substrates (EEF1A Y254, RPS27 S11, RPL27 S34, etc.), which are involved in translation and ribosome biogenesis were upregulated in proteomic subtype GP3 (Fig. 4G). Furthermore, we observed that the proteomic subtype GP3 has the highest ribosome biogenesis pathway score (Fig. 4H, Kruskal–Wallis test, p < 0.05) through sample-specific gene set enrichment analysis (ssGSEA, Methods). Elevated ribosome biogenesis pathway signaling was associated with decreased survival (Fig. 4I and 4J, log rank test, p < 0.05), either in all proteomic subtypes or specifically in the GP1 and GP2 subtypes. The association of the ribosome biogenesis pathway with worse prognostic outcomes was also confirmed in the Chicago T1BC cohort (Figure S4E, log rank test, p < 0.05). We then classified GP1/2 subtypes into GP1/2 with a high ribosome biogenesis pathway score (GP1/2-RH) and GP1/2 with a low ribosome biogenesis pathway score (GP1/2-RL), based on the survival group presented in Fig. 4J. Integrative analysis of the BCG-NR and BCG-R groups, three proteomic subtypes, and the expression level of TMB and TLR3, revealed a high consistency between the GP1-RL and GP2-RL proteomic subtypes ,the BCG-R group and the groups with a high level of TMB or TLR3 (Fig. 4K; Table S4). These results suggested that patients with a low level of TMB and TLR3 but a high level of ribosome biogenesis might not benefit form BCG treatment and could potentially benefit from inhibitors of ribosome biogenesis.

Construction and validation of the predictive models for BCG therapy

Having described the genomic, proteomic, and phosphoproteomic features associated with response to BCG therapy, we next used a machine learning framework to integrate these features and create a predictive model of the BCG response (Fig. 5A). We applied Wilcoxon rank-sum tests with a Benjamini-Hochberg adjusted p value cutoff (FDR < 0.05) and identified 125 differentially expressed proteins between BCG-NR tumors and BCG-R tumors (Figure S5A, Wilcoxon rank-sum test, FDR < 0.05; Table S5). We utilized stepwise logistic regression, a method that is resistant to noise and overfitting, to identify a subset of proteins (BCG-R/BCG-NR-sig) that differentiate between BCG-NR and BCG-R. To train and subsequently test the model, the samples were partitioned based on their type (i.e., R or NR). The training set consisted of 75% of the samples, while the remaining 25% were used as the testing set. Based on BCG-R/BCG-NR-sig (N = 6, Figure S5B; Table S5), we applied 10-fold cross-validation on the training set, which resulted in a predictive model with a mean receiver operating characteristic-area under the curve (ROC-AUC) of 0.819 (Fig. 5B). When applied to the testing set samples, the predictive model achieved a ROC-AUC of 0.833, a specificity of 0.85, and a sensitivity of 0.83 (Fig. 5B and Figure S5C).

To assess the accuracy of the predictive model for BCG response, we employed the targeted MS approach, parallel reaction monitoring (PRM) assays, which has been adopted in classifier’s validation in recent proteomic research [43, 44], to quantify the proteins of the classifier in tumor tissues from an independent validation cohort consisting of 50 NMIBC patients who underwent BCG therapy (BCG-R, N = 38; BCG-NR, N = 12). We selected a set of target peptides unique to these proteins (TLR3, PSMA2, OS9, LSM3, PHAX, and RBCK1) based on the library search results (Table S5). Through the PRM quantification, we observed a significant differential expression of these proteins between BCG-NR and BCG-R (Fig. 5C, Wilcoxon rank-sum test, FDR < 0.05), with a mean ROC-AUC of 0.854 (Fig. 5C). Collectively, the proteomic classifier shows promise as a potential predictive model for distinguishing between BCG-NR and BCG-R.

To further improve the predicted performance of the classifier, we enrolled the TMB, a genomic feature associated with the response to BCG therapy, in addition to the proteome features, resulting in a multi-omics classifier. The multi-omics classifier had a mean ROC-AUC of 0.909 and 0.867 on the 75% training set and 25% testing set, respectively, outperforming compared to the proteomic classifier in the discovery cohort (Fig. 5D and Figure S5D). The sensitivity and specificity of the testing set were both 1, indicating a significant improvement over the proteomic classifier (Fig. 5D). In summary, we used an ensemble approach that input multi-omics features to identify predictors for the BCG Responders. This predictive model offers an opportunity to expedite the translation of basic research into more precise diagnosis and treatment in the clinic.

Progression Clock of Patients with NMIBC received BCG therapy

Survival analysis showed that the BCG-R group had a higher survival rate than the BCG-NR group (Fig. 6A-B, log rank test, p < 0.05). When making decisions about treatment and follow-up, it is crucial to consider disease progression as the primary factor. Predicting the progression of NMIBC after BCG therapy is a challenging task in the medical field, but it plays a significant role in planning disease treatment. To predict the short- and long-term probabilities of disease recurrence and progression, we sought to build time-dependent models based on proteome data. Firstly, to investigate the association between protein expression and BCG therapy, we integrated the 73 pre-BCG samples and performed weighted gene correlation network analysis (WGCNA) (Figure S6A). The separation of proteomic profiles into eigengene modules identified a specific group of modules that positively associated with PFS, papillary carcinoma, and risk score (Fig. 6C). Among these modules, MEblack was significantly correlated with PFS time and OS time (p < 0.05). The genes within this module were enriched in apoptosis and immune-related pathways, such as the regulation of immune effector process (Fig. 6D; S6B).

To further explore the prognostic value of proteins, we constructed a univariate Cox hazard model to filter the proteins in the MEblack module. The result showed that 10 proteins, including ATP6A2, PUS1, HAUS3, GOLM1, SDF2L1, TBCD, LMNB2, DNAJC3, OSTF1, and TXNDC5, met the predefined criteria (Fig. 6E). Subsequently, we established a PFS score using a multivariate Cox proportional hazard regression model that incorporated these 10 proteins (Fig. 6F, Methods).

To further assess the potential clinical utility of the model, we conducted a test to determine if the predicting protein model could be used to predict the progression time in patients with NMIBC. The progression clock operating characteristic analysis of the PFS score showed that the AUC was 0.902, 0.832, and 0.803 at 1-, 3-, and 5-year, respectively (Fig. 6G). Subsequently, we integrated the ten proteins to construct a nomogram (Figure S6C). The calibrated plots demonstrated that the predicted 1-, 3-, and 5-year survival probabilities of the nomogram performed well in our NMIBC cohort (Figure S6D). Additionally, decision curve analysis was used to assess the predicting the advantages of nomogram (Figure S6E), and decision curve analysis (DCA) displayed a net benefit of the nomogram. We further evaluated the performance of our model in another independent validation cohort (N = 74), where the AUC was 0.95, 0.81, and 0.87 at 1-, 3-, and 5-year, respectively (Figure S6F). These results indicated that the strength and the high degree discriminative power of our model in predicting the prognosis of NMIBC patients treated with BCG.

TLR3-Mediated Activation of NF-κB Pathway and Immune Checkpoint Modulation in Bladder Cancer: Insights from Experimental Analyses

Consistent with previous research findings [45, 46], our correlation analyses revealed a significant positive association between NF-κB1 and TLR3 (Figure S3G). To further investigate the involvement of the NF-κB signaling pathway in TLR3-mediated regulation in bladder cancer, we utilized poly I:C (a TLR3 activator) and TLR3 overexpression cell lines. Immunofluorescence analysis revealed co-localization of TLR3 and NF-κB1 molecules, suggesting a potential interaction between them (Fig. 7A). When TLR3 was activated by poly(I:C) or overexpressed, there was an upregulation of phosphorylation levels of p65, indicating activation of the NF-κB pathway (Fig. 7B). This suggested that the NF-κB pathway was activated in response to TLR3 stimulation by poly(I:C) or overexpression. It was crucial to highlight that these findings not only shed light on the broader regulatory role of TLR3 but also have direct implications for BCG therapy sensitivity in BC. The observed activation of the NF-κB pathway, a key player in inflammatory responses and cellular processes [46, 47], suggested that TLR3 may contribute to modulating the immune microenvironment and, consequently, impact the efficacy of BCG treatment. This insight strengthens the rationale for considering TLR3 as a potential target for enhancing the sensitivity of bladder cancer to BCG therapy. Furthermore, our correlation analysis also unveiled a positive association between TLR3 expression and various MHC molecules (Figure S3D). Following TLR3 overexpression, there was indeed an upregulation of HLA-I and PD-L1 in tumor cells (Fig. 7C; Fig. 7D; Fig. 7E). The upregulation of HLA-1 is crucial for antigen presentation, and the concomitant increase in PD-L1 may have implications for immunotherapy strategies [48, 49]. In conclusion, TLR3 may play a role in modulating BCG responses by regulating the antigen presentation process.

TLR3-Mediated Regulation of Cytokine Secretion in Bladder Cancer Cell Lines

To gain a deeper understanding of the functionality and regulatory mechanisms of TLR3, we examined alterations in the cytokine profile. We found that TLR3 overexpression led to increased secretion of IL-4 and IL-10 in both T24 and 5637 cell lines (Fig. 7G). However, it was noteworthy that RT4 cells did not exhibit significant alterations in cytokine secretion when TLR3 was overexpressed. These findings suggested that the response to TLR3 modulation is specific to the type of cell. The heightened secretion of IL-4 and IL-10 in T24 and 5637 cells could potentially contribute to immunomodulation within the tumor microenvironment [50], fostering an environment conducive to tumor progression. To sum up, these results emphasize the importance of considering cell-specific responses when exploring the immunomodulatory effects of TLR3 overexpression in bladder cancer, offering valuable insights for the development of targeted therapeutic strategies tailored to specific cellular contexts.

As a pattern recognition receptor (PRR), TLR3 plays a pivotal role in innate immunity [51]. Therefore, our study aimed to investigate the role of TLR3 in the activation of innate immunity following BCG stimulation. A previous study showed that anthracycline chemotherapy drugs activate TLR3, stimulating the production of type I interferon (IFN) by breast cancer tumor cells [52]. This self-secreted IFN, upon binding to IFN receptors on tumor cells, exerted autocrine-mediated anti-tumor effects. Our results revealed that TLR3 overexpression increases the secretion of IFN-β (Fig. 7F) and enhances cellular apoptosis after BCG treatment (Figure S7C). Importantly, when IFN receptors were blocked, the observed increase in apoptosis was attenuated (Figure S7C). These findings provided a partial explanation for the heightened sensitivity to BCG therapy observed with TLR3 overexpression. Furthermore, the heightened secretion of IFN-β, along with an augmented apoptotic response to BCG, emphasizes the potential influence of TLR3 modulation on the cellular dynamics of bladder cancer. This suggested a promising historical therapeutic approach for manipulating the immune response against tumor cells [53].

Enhanced Sensitivity of TLR-High Expressing Organoids to BCG

To comprehensively elucidate the impact of TLR3 on BCG treatment sensitivity, we explored the association between TLR3 expression and BCG responsiveness in a BC organoid model. Utilizing tumor tissues from 10 bladder cancer patients, we established 10 patient-derived organoids (PDOs), comprising 5 PDOs exhibiting high TLR3 expression and 5 with low expression (Figure S7B). Following a 7-day BCG treatment, we observed that organoids with high TLR3 expression exhibited increased sensitivity to BCG (Fig. 7H; Figure S7D). Morphologically, the TLR-high-expressing organoids displayed significant shrinkage after BCG treatment, along with a more pronounced decrease in ATP activity (Fig. 7H; Figure S7D). In summary, our investigation into the association between TLR3 expression and BCG responsiveness in a BC organoid model provides intriguing insights. The enhanced sensitivity of TLR-high-expressing organoids to past BCG treatments highlights the historical importance of TLR signaling in the response to immunotherapeutic interventions. Modulating TLR expression or activity may enhance the efficacy of BCG therapy, and further investigations into these pathways may unveil novel targets for therapeutic intervention in BC.

We performed an integrative multi-omics analysis of NMIBC tumors from 86 patients treated with BCG. Our analysis provides a more comprehensive understanding of the molecular characteristics of BCG responses and may contribute to explaining BCG-resistance.

WES analysis showed that well-known BC genes, such as KMT2D, KDM6A, TP53 and FGFR3, were not associated with response to BCG immunotherapy or improved PFS. Mutational signature analysis indicated that APOBEC, which had the highest TMB, was associated with the best prognosis. It has been reported that TMB, a potential indicator of tumor immunogenicity, is an emerging predictive biomarker of response to ICI treatments. Many studies have demonstrated a strong correlation between TMB levels and clinical outcomes [54–58]. In this study, we investigated the predictive value of TMB in Chinese patients with NMIBC following BCG therapy and identified additional biomarkers that might help refine the association between TMB and response. Furthermore, the limited understanding of the mechanism of BCG action does not fully explain the correlation between TMB and BCG. We found that high TMB was correlated with high STAT1 activity, which might activate antigen processing and presentation through PSMB9. Many studies have reported that STAT1 can upregulate the constitutive transcription of several genes, mainly involved in the presentation of tumor antigens (MHC class, CIITA, and LMP2) [59, 60]. Increased expression of these genes leads to more efficient presentation of antigens in the context of class I MHC molecules [61]. Additionally, PSMB9 is involved in antigen processing to generate class I binding peptides, and the MHC class I restricted pathway of antigen processing allows the presentation of intracellular antigens to cytotoxic T lymphocytes [62–64]. Higher expression levels of PSMB9 might enhance the therapeutic effect of BCG, which should be validated in future studies. Taken together, our data provided new insights into the relationship between high TMB and BCG response, serving as a resource for further mechanistic and functional research related to BCG therapy.

Genetic changes are always associated with bladder cancer progression. Our data also revealed significant correlations between arm events, such as Amp 1q21.1, Amp1q23.3, and Del 4q35.1, and clinical outcomes. It is assumed that circumscribed high-level amplification/deletion is linked to the expression of oncogenes. We found that higher expression of TLR3 (4q35.1) appeared to have better prognostic outcomes, and TLRs have been reported to be markedly decreased in BC samples than normal urothelial samples [33]. Recently, TLR agonists have been considered as the better adjuvants for vaccine immunotherapy in cancer therapy, and the TLR3 agonist, poly(I:C) has shown potent direct cytotoxic and immunomodulating activity on low-grade human BC cell lines [65]. Additionally, the integration of TMB and TLR3 expression to identify potential responders to BCG therapy remains unclear. Combined analysis showed that patients with high TMB and high TLR3 expression had the best prognosis, indicating that TMB can further stratify the outcomes of BCG therapy for NMIBC patients within each clinically relevant TLR3 expression group. Thus, our data suggested that high TMB levels were associated with improved BCG efficacy across different TLR3 expression subgroups, which holds significant implications. Detecting TLR3 and TMB through WES may help distinguish patients who will respond effectively to BCG therapy.

In this study, distinct proteomic subtypes were identified, each with unique molecular features and clinical outcomes, and associated with BCG therapy. The patients in the GP1 proteomic subtype had the most favorable prognosis and were more likely to belong to the LumGU subtype, previously identified by Hurst et al. in the stage T1BC Chicago cohort [18]. The GP1 proteomic subtype exhibited increased expression of genes related to phagosome and antigen processing and presentation, which is vital for MHC cytotoxic T cell killing of cancer cells [66, 67]. In contrast, GP3 proteomic subtype had the poorest prognosis and the highest ribosome biogenesis pathway score. Recent studies have indicated that dysregulation of ribosome biogenesis was associated with malignant tumor behavior, therapeutic resistance, and clinical outcomes [68, 69]. Interestingly, we found that high expression of ribosome biogenesis was significantly associated with a poor prognosis, even in the GP1 and GP2 subtypes with favorable prognosis (GP1/2-RH). In addition, we observed that the GP1-RH and GP2-RH proteomic subtypes were highly consistent with the BCG-NR group as well as the groups with low levels of TMB or TLR3. These results suggested that patients with low levels of TMB and TLR3 but high levels of ribosome biogenesis may not benefit from BCG treatment but could potentially benefit from inhibitors targeting ribosome biogenesis. Thus, targeting the ribosome could be a viable therapeutic option for patients who have not responded to BCG therapy.

For over 40 years, BCG has been the most effective intravesical therapy for high-risk disease. yet pretreatment biomarkers that can reliably predict which patients will benefit from BCG treatment are still needed [70]. Clinically localized NMIBC patients deserve the most accurate and tailored predictions available, which current staging systems do not provide. As previously reported models for predicting the response to BCG were almost constructed based on transcriptome data [71], constructing a classifier to predict the drug response based on multi-omics data provides an alternative perspective. In this study, we first employed proteomic technology implemented with machine learning statistics to develop predictive models for BCG response and non-response. The identified signature proteins were validated using PRM quantification in an independent cohort, effectively distinguishing between BCG-NR and BCG-R. Moreover, when we incorporated the genomic feature TMB related to the response to BCG treatment, a multi-omics classifier was generated, which outperformed the proteomic classifier. In conclusion, this classifier can identify NMIBC patients who failed BCG treatment, regardless of traditional clinical parameters such as tumor size, tumor stage or the number of tumors. This may allow individual risk stratification and potentially may influence the management of NMIBC patients.

BCG immunotherapy markedly improves the outcome of high-risk NMIBC [9]. Nevertheless, there is a subgroup of patients who fail to response to the treatment and are at risk for disease progression. In these patients, tools to predict the time of progression in NMIBC patients after BCG therapy would be significant. Despite extensive research currently, molecular markers for recurrence and progression have not yet been widely used in clinical practice [16], and there is still a lack of the molecular markers for predicting the time of progression. In this study, we developed a model based on progression-related proteins to predict the time of BCG failure using proteomic data. We also further proved that this model is a reliable indicator for predicting 1-, 3-, and 5-year survival of NMIBC patients treated with BCG therapy, and its effectiveness was also verified in an independent cohort. These results indicated that the features of progression-related protein could be a potential predictive model to predict the progression time on the basis of tissue proteome from NMIBC patients treated with BCG therapy. This could lay the foundation for further improvements towards a personalized approach to BCG therapy in NMIBC and may have value in clinical applications.

In conclusion, we reported an integrative multi-omics analysis of NMIBC tumors from 86 NMIBC patients treated with BCG. This is the first study to provide a comprehensive multi-omics landscape of NMIBC treated with BCG in an Asian population. Through proteomic classification, we identified three molecular subtypes in pre-BCG NMIBC that have prognostic relevance, which provides a framework for biomarker discovery to optimize the current clinical management of patients with NMIBC treated with BCG therapy. More importantly, current evidence supported the use of TMB high and TLR3 high as biomarkers for BCG treatment in NMIBC. However, future studies should focus on assessment of TMB from target sequencing and the investigations of specific activity of BCG therapy in TMB high and TLR3 high tumors before these biomarkers can be widely implemented in clinical practice.

Patient cohort

In this study, patients with treatment- naïve histologically confirmed NMIBC treated with TURBT and induction intravesical BCG between 2012 and 2021 were identified from the Department of Urology of Fudan University Shanghai Cancer Center (FUSCC, Shanghai, China). Studies were performed in accordance with the Declaration of Helsinki. All patients with T1 disease underwent a restaging TURBT. Patients were selected based on the availability of bladder tissue specimen pre-BCG treatment with enough tumor available to obtain a core. Specially, pre-BCG treatment samples were obtained from primary, incident tumors 1–3 months prior to starting a 6-week induction course of BCG. Cystoscopic evaluations were performed 3 months after the first BCG instillation, and those responding to treatment were then evaluated with cystoscopy and urine cytology every 3 months for 2 years, 6 months to 5 years, and annually up to 10 years at the discretion of the treating urologist. All pathologic evaluations were performed at FUSCC by a genitourinary pathologist. BCG response was defined according to standard definitions: (1) patients with a minimum follow-up of 2 years after diagnosis were defined as “BCG Responders (BCG-R)”. (2) “BCG non-responders (BCG-NR)” included BCG unresponsive, BCG relapsing, and BCG progressor. BCG unresponsive patients were those who had persistent high-grade T1 disease at the initial 3-month cystoscopy, or had relapsed high-grade NMIBC with or without CIS within 6 months of the last exposure to BCG. BCG relapsing patients had recurrent high-risk (high-grade) NMIBC after a prior complete response and did not fulfill the BCG-unresponsive definition [72, 73]. Patients who received an inadequate BCG induction (less than five of six courses) and those with incomplete follow-up information were excluded.

Sample preparation

Tumor tissue samples were collected from the FUSCC tissue bank immediately after surgery. The samples were collected within 30 minutes after resection, and were immediately transferred into sterile freezing vials and snap frozen in liquid nitrogen. The tumor tissues were then cut into approximately 0.5 cm³ pieces under a temperature of -40°C and stored at -80°C until further use. Histologic sections were obtained from the top and bottom portions of the tumor tissues and stained with Hematoxylin and eosin (H&E) for review. In addition, precancerous lesions and divergent histological variant tumors of patients were scraped according to H&E staining. Tumor samples were required to contain an average of 70% tumor cell nuclei with necrosis equal to or less than 20% for inclusion in the study. Each sample was assigned a new research ID, and the patient’s name or medical record number used during hospitalization was de-identified.

Pathology Review

All samples were systematically evaluated to confirm the histopathologic diagnosis and any variant histology according to the World Health Organization (WHO) classification by three expert genitourinary pathologists. Additionally, the tumor samples were also assessed for tumor content, presence and extent of tumor necrosis, and signs of invasion into the lamina propria or muscularis propria. Tumor samples were also evaluated for the presence and extent of inflammatory infiltrates, as well as for the type of the infiltrating cells (lymphocytes, neutrophils, eosinophils, histiocytes, plasma cells) in the tumor microenvironment. Any non-concordant diagnoses among the three pathologists were re-reviewed, and a resolution was reached following discussion. All the information is included in Supplementary Data 1.

Whole-exome Sequencing

DNA extraction

For WES, DNA was extracted from fresh or frozen tumor tissue samples, while matched germline DNA was obtained from corresponding blood samples. DNA from both the tumor tissues samples and blood samples was extracted according to the manufacturer’s instructions of a DNeasy Blood & Tissue Kit (QIAGEN, Hilden, Germany). The quality of the isolated and contaminated samples was verified using the following methods: (i) monitoring DNA degradation and contamination using 1% agarose gels; and (ii) measuring DNA concentration with the Qubit® DNA Assay Kit in a Qubit® 2.0 Fluorimeter (Invitrogen, CA, USA).

Library preparation

An input of 0.6 µg genomic DNA per sample was used for the DNA preparation. The DNA sequencing libraries were generated following the manufacturer's recommendations using an Agilent SureSelect Human All Exon kit (Agilent Technologies, CA, USA) following the manufacturer’s recommendations, following which index codes were added to each sample. Briefly, fragmentation was carried out using a Covaris M220 system to generate 200–300 bp fragments. The remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. Following adenylation of the 3’ ends of DNA fragments, adapter oligonucleotides were ligated. DNA fragments with adapter molecules ligated at both ends were selectively enriched using a PCR reaction. Following the PCR reaction, the libraries were then hybridized with t a biotin-labeled probe in the liquid phase and captured using magnetic beads with streptomycin to capture the exons of genes. Enriched libraries were subjected to a PCR reaction to add index tags for sequencing preparation. The products were purified using an AMPure XP system (Beckman Coulter, Beverly, USA) and quantified using the Agilent high sensitivity DNA assay on the Agilent Bioanalyzer 2100 system.

Clustering and Sequencing

Clustering of index-coded samples was performed using a cBot Cluster Generation System with a Hiseq PE Cluster Kit (Illumina) according to the manufacturer’s instructions. After cluster generation, the DNA libraries were sequenced on the Illumina HiSeq platform and 150 bp paired-end reads were generated.

WES data analysis

Quality Control

The original fluorescence image files obtained from the Hiseq platform were transformed to short reads (raw data) through base calling. These short reads were then recorded in FASTQ format, which contains both sequence information and corresponding sequencing quality information. To ensure reliable bioinformatics analysis, it is crucial to address sequence artifacts, including reads containing adapter contamination, low-quality nucleotides, and unrecognizable nucleotides (N). Hence, quality control is an essential step that must be applied to guarantee meaningful downstream analysis.

The data processing steps were as follows:

Paired reads were discarded if either read contained adapter contamination (> 10 nucleotides aligned to the adapter, allowing ≤ 10% mismatches).
Paired reads were discarded if more than 10% of bases are uncertain.
Paired reads were discarded if the proportion of low-quality (Phred quality < 5) bases is either read was over 50%.

All downstream bioinformatics analyses were based on high-quality clean data, which were retained after these steps. At the same time, QC statistics including total read number, raw data, raw depth, sequencing error rate, percentage of reads with Q30 (the percentage of bases with Phred-scaled quality scores greater than 30), and GC content distribution were calculated and summarized.

Reads Mapping to Reference Sequence

Valid sequencing data were mapped to the reference human genome (UCSC Genome Brower, hg38) using Burrows-Wheeler Aligner (BWA) software [74] to obtain the original mapping results stored in BAM format. If one read, or one paired read, were mapped to multiple positions, the strategy adopted by the BWA was to choose the most likely placement. If two or more most likely placements were present, the BWA picked one randomly. Then, SAMtools [75] and Picard (http://broadinstitute.github.io/picard/) were used to sort BAM files and perform duplicate marking, local realignment, and base quality recalibration to generate final BAM files for computation of the sequence coverage and depth. The mapping step was very difficult due to mismatches, including true mutations and sequencing errors, and duplicates resulting from PCR amplification. These duplicate reads were uninformative and should not be considered as evidence for variants. We used Picard to mark these duplicates for the follow-up analysis.

Variant Calling

Somatic variants were then called, utilizing VarScan v2.3.8 [76] MuTect v1.1.7 [77], and InVEX (http://www.broadinstitute.org/software/invex/). The following filters were applied to get variant cells of high confidence: Remove mutations with coverage less than 10×; Remove variant sites in dbSNP and with mutant allele frequency (MAF) > 0.05 in the 1000 Genomes databases (1000 Genomes Project Consortium) and the Novo-Zhonghua (in-house unrelated healthy individual database), but include sites with MAF ≥ 0.05 with COSMIC evidence (http://cancer.sanger.ac.uk/cosmic) [78, 79]; All variants must be called by 2 or more callers; All variations must be exonic; Retain the nonsynonymous SNVs if the functional predictions by PolyPhen-2, SIFT, MutationTaster and CADD all show the SNV is not benign [80, 81]; Retain genes identified by Cancer Gene Census (CGC, http://www.sanger.ac.uk/science/data/ cancer-gene-census). The genes with mutations were analyzed by Fisher's exact test to compare the two proteomic groups.

Copy-Number Analysis

CNAs were identified using GATK’s (GATK 4) Best Practice somatic CNA calling pipeline. The results of this pipeline, segment files of every 1000 were then input into GISTIC2 [82], to identify significantly amplified or deleted focal-level and arm-level events, with a q-value < 0.1 considered significant. A log2 ratio cut-off 1 was used to define SCNA amplification and deletion. We further summarize the arm-level copy number change based on a weighted sum approach [83], in which the segment-level log2 copy ratios for all segments in the given arm were added up with the length of each segment being weighted. To exclude false positives as much as possible, relatively stringent cut-off thresholds were used with the following parameters: -ta 0.1 -tb 0.1 -brlen 0.98 -conf 0.9. Other parameters were the same as their default values.

Co-occurrence and Mutual Exclusivity Analysis of Mutations

Co-occurrence and mutually exclusive mutated genes were detected using Fisher’s exact test to determine the co-occurrence and mutually exclusively of significantly mutated genes in our mutational dataset.

Analysis of Significantly Mutated Genes

Filtered mutations (including SNV and indel) were further used to identify significantly mutated genes using MutSigCV (https://software.broadinstitute.org/cancer/cga/mutsig, version 1.4) with default parameters. The final MutSigCV P values were then converted to q-values using the method of Benjamini and Hochberg [84]. Genes with q-value < = 0.1 were declared to be significantly mutated.

Mutational signature analysis using the sigminer approach

Mutational signatures were jointly inferred for 92 tumors using the R package sigminer [85]. The sigminer approach (https://github.com/ShixiangWang/sigminer) was used to extract the underlying mutational signatures. The input data consisted of 96 mutation vectors (or contexts) generated by somatic SNVs based on six base substitutions (C > A, C > G, C > T, T > A, T > C, and T > G) within 16 possible combinations of neighboring bases for each substitution. These mutation vectors were used to determine their contributions to the observed mutations. Sigminer using a non-negative matrix factorization (NMF) approach was applied to decipher the 96 × 92 (i.e., mutational context-by-sample) matrix for the 30 known COSMIC cancer signatures (https://cancer.sanger.ac.uk/cosmic/signatures) and infer their exposure contributions.

Mutational signature analysis using the deconstructSigs approach

The mutational signature of each sample was deconstructed using the deconstructSigs approach [86] and its R package (deconstructSigs v1.8.0) with default parameters. Thirty COSMIC cancer signatures were considered and their contributions (weights) in each patient were normalized between 0 and 1, and signatures with a weight below 0.08 were filtered out.

Functional Annotation

Functional annotation plays a vital role in elucidating the link between genetic variations and diseases. ANNOVAR was performed to annotate the variant call format (VCF) obtained in a previous study [87]. The dbSNP, 1000 Genome, and other related databases were used to characterize the detected variants. Given the significance of exonic variants, gene transcript annotation databases, such as Consensus CDS, RefSeq, Ensembl, and UCSC, were also included in the determination of amino acid alterations. Annotation content contained the variant position, variant type, and conservative prediction, and more. These annotation results would help locate disease-causal mutants. The details of the annotation are provided in the supplementary material.

TMB

TMB was defined as the number of somatic mutations (including base substitutions and indels) in the coding region. Synonymous alterations were also counted [88]. To calculate the TMB, the total number of mutations counted was divided by the size of the coding sequence region of the xGen Exosome Research Panel (Integrated DNA Technologies, IDT). The WES target region was 33.52M. In this study, a threshold value of10 was adopted to differentiate between TMB high (TMB-H) and TMB low (TMB-L) cases.

Proteomic and phosphoproteomic analysis

Protein extraction and trypsin digestion

Frozen samples were collected and washed three times with phosphate buffer saline (PBS) to remove blood and debris. The samples were then minced and treated with lysis buffer (8M urea, 100 mM Tris-HCl, pH 8.0) containing protease and phosphatase inhibitors (Thermo Scientific, USA) and then sonicated for 1 minute (3s on and 3s off, amplitude 25%). The resulting lysates were centrifuged at 14,000g for 10 minutes and the supernatants were collected as whole-tissue extracts. Protein concentrations were determined by the Bradford protein assay (TaKaRa, T9310A). Extracts (100 µg protein) were reduced with 10 mM dithiothreitol at 56°C for 30 minutes and alkylated with 10 mM iodoacetamide at room temperature in the dark for 30 minutes. Finally, the samples were digested with trypsin using a filter-aided sample preparation method [89].

The enrichment of phosphorylated peptides

For the phosphoproteomic analysis, peptides were extracted from the FFPE slides following the same protocol as previously described for “Protein extraction and trypsin digestion”. Tryptic peptides were then used for phosphopeptide enrichment using a High-Select Fe-NTA kit (Thermo Fisher Scientific, Rockford, IL, USA, #A32992) according to the kit manual and a previous report [90] with some modifications. In brief, the peptides were suspended in the binding/washing buffer (contained in the enrichment kit) and mixed with the equilibrated resins. The peptide-resin mixture was incubated for 30 minutes with three gentle blows at room temperature. Following incubation, the resins were washed three times with the binding/washing buffer and twice with water. The enriched peptides were eluted with the elution buffer (contained in the enrichment kit) and immediately dried using a speed-vac at 45°C for mass spectrometry analysis.

LC-MS/MS analysis

For the analysis of proteome profiling samples and phosphoproteomic samples, peptides were analyzed on an Orbitrap Fusion Lumos Hybrid Quadrupole-Orbitrap Mass Spectrometer (Thermo Fisher Scientific) coupled with a high-performance liquid chromatography system (EASY nLC 1200, Thermo Fisher Scientific). Dried peptide samples were re-dissolved in Solvent A (0.1% formic acid in water) and loaded onto a 2-cm self-packed trap column (100 µm inner diameter, 3 µm ReproSil-Pur C18-AQ beads, Dr Maisch GmbH) using Solvent A. The peptides were then separated on a 150 µm-inner-diameter column with a length of 30 cm (1.9 µm ReproSil-Pur C18-AQ beads, Dr Maisch GmbH) over a 150 mins gradient (Solvent A: 0.1% Formic acid in water; Solvent B: 0.1% Formic acid in 80% ACN) at a constant flow rate of 600 nL/min (0-150 mins, 0min, 4% B; 0–10 mins, 4%-15% B; 10–125 mins, 15%-30% B; 125–140 mins, 30%-50% B; 140–141 mins, 50%-100% B; 141–150 mins, 100% B). The eluted peptides were ionized at 2 kV and introduced into the mass spectrometer. Mass spectrometry was performed in data-dependent acquisition mode. For the MS1 Spectra full scan, ions with m/z ranging from 300 to 1,400 were acquired by an Orbitrap mass analyzer at a high resolution of 120,000. The automatic gain control (AGC) target value was set to 3E + 06 with a maximal ion injection time of 80 ms. MS2 spectral acquisition was performed in the ion trap in a rapid speed mode with a cycletime of 1.5s. Precursor ions were selected and fragmented with higher energy collision dissociation (HCD) with a normalized collision energy of 30%. Fragment ions were analyzed by an ion trap mass analyzer with an AGC target set at 5E + 04 and the maximal ion injection time of MS2 was 20 ms. Peptides that triggered MS/MS scans were dynamically excluded from further MS/MS scans for 18 s. The coefficient of variation (CV) values on FAIMS were − 45V and − 65V.

MS database searching

Peptide and protein identification

MS raw files were processed using the “Firmiana” (a one-stop proteomic cloud platform) [91] against the human National Center for Biotechnology Information (NCBI) RefSeq protein database, with the Mascot 2.4 software (Matrix Science Inc., London, UK). The search parameters included a maximum of 2 missed cleavages, mass tolerances of 20 ppm for precursor ions and 0.05 ppm for production ions and fixed modification of cysteine carbamidomethylation. Variable modifications considered were N-acetylation and methionine oxidation. To ensure the quality control of protein identification, a target-decoy-based strategy was implemented to maintain a FDR of less than 1% for both peptides and proteins. The program percolator was used to obtain the probability value (q-value) and showed that the FDR (measured by the decoy hits) of every peptide-spectrum match (PSM) was lower than 1%. Peptides shorter than seven amino acids were excluded, and a cutoff ion score of 20 was set for peptide identification. To ensure stringent quality control, PSMs from all fractions were combined. The q-values of both target and decoy peptide sequences were dynamically increased employing the parsimony principle until the corresponding protein FDR was less than 1%. Finally, to reduce the false positive rate, proteins with at least two unique peptides were selected for further investigation.

Label-free-based MS Quantification of Proteins

For the calculation of protein abundance, the non-redundant peptide list was used to assemble proteins following the parsimony principle. Protein abundance was then estimated by a traditional label-free, intensity-based absolute quantification (iBAQ) algorithm, which divided the protein abundance (derived from identified peptide intensities) by the number of theoretically observable peptides. Match between runs [92] was used to improve parallelism across different tissues from 86 patients. We built a dynamic regression function based on commonly identified peptides in different tissues. According to the correlation value, R2, Firmiana chooses a linear or quadratic function for regression to calculate the retention time (RT) of the corresponding hidden peptides, and check the existence of the extracted ion chromatogram (XIC) based on the m/z and calculated RT. The program evaluated the peak area values of the existing XICs which were calculated as parts of the corresponding proteins. Proteins with at least two unique peptides with a 1% FDR at the peptide level were selected for further analysis. Then, the fraction of total [48], a relative quantification value that was defined as a protein’s iBAQ divided by the total iBAQ of all identified proteins in one experiment, was calculated as the normalized abundance of a particular protein among experiments. Finally, the fraction of total [48] was further multiplied by 10⁵ for ease of presentation and any missing values were assigned as 10^− 5 [93].

Batch effect analysis

Hierarchical clustering, dip statistic test and principal component analyses were implemented in R v.3.5.1 to assess batch effects in our proteome dataset with respect to the following two variables: batch identity and sample type (pre-BCG and normal samples). For the hierarchical clustering analysis, pairwise Spearman’s correlation coefficients of samples that passed quality control were investigated. Samples of the same type exhibited high similarity, whereas samples of different types exhibited clear differences. There was no evident association between batch identity and correlation coefficients. The density plot of the normalized intensities of the identified proteins in each sample showed that all samples passed quality control as expected from a unimodal distribution (dip statistic test). The results of principal component analysis showed that batch effects were negligible for batch identity but significant for the sample types.

Quality control of the MS data

For quality control of MS, the HEK293T cell (National Infrastructure Cell Line Resource) lysates were measured every three days to set the quality-control standard. The quality-control standard was then digested and analyzed using the same method and conditions as the 12 samples. Pairwise Spearman’s correlation coefficient was calculated for all quality-control runs using the R package corrplot (v0.84), and the results are shown in Figure S1F. The average correlation coefficient among the standards was 0.90, while the maximum and minimum values were 0.95 and 0.86, respectively.

Targeted PRM analysis

To evaluate the accuracy of the classifiers for predicting BCG therapy, we designed a PRM strategy to quantify the classifier proteins in plasma samples from an independent cohort composed of 50 NMIBC patients treated with BCG. From the library search results, a set of target peptides that were unique to the classifier proteins were selected. Besides, house-keeping proteins such as VCP, RPLP0, and PSMB4, were also included as the reference. Equal amounts of tissue from each sample in the validate cohort with 50 NMIBC patients were digested as described in the profiling preparation section. The peptide samples were then injected into the Q Exactive HF-X Hybrid Quadrupole-Orbitrap Mass Spectrometer (Thermo Scientific) operating in PRM mode with quadrupole isolation and HCD fragmentation. The mobile phase buffers used were the same as those mentioned in the DIA method section. Full MS mode was measured at a resolution of 60,000 with an AGC target value of 3E6 and maximum IT of 20ms, scanning the range of 300–1400 m/z. Target precursors were then isolated through an m/z window of 1.2 Th, followed by fragmentation at 27% normalized collision energy. The product ions were scanned with a resolution of 15,000, an AGC target value of 1E6 or a maximum injection time of 25ms. The raw data was searched by Skyline-daily (4.2.1.19004, University of Washington, USA). The proteins were quantified based on the fragment total area reported by Skyline-daily. We selected peptides and tested their signal stability and peak shape in the pool sample for final quantification, following the ranking offered by Skyline.

Quantification and Statistical Analysis

Missing Value Imputation

For the proteomic and phosphoproteomic data, FOTs were multiplied by 1E5 for quantification. Proteins and phosphoproteins detected in more than 30% of the samples underwent missing value imputation. The missing values were imputed with 1E-5 and, if necessary, log2 transformed. This method has also been previously applied in other published proteomic and phosphoproteomic studies [93, 94].

Differential protein analysis

Proteins expressed in more than 30% of the samples were selected for differential expression analysis. The Wilcoxon rank-sum test was used to examine whether proteins were differentially expressed between BCG-R (n = 52 samples) and BCG-NR (n = 21 samples), or patients with different mutation statuses and CNA statuses. Upregulated or downregulated proteins are defined as proteins that show differential expression in one group compared with the other group (Wilcoxon rank-sum test, FDR < 0.05, Fold change > 2 or < 0.5). The Kruskal-Wallis test was used to test whether genes were differentially expressed among different tissue types or other subgroups. To account for multiple-testing, the P values were adjusted using the Benjamini-Hochberg FDR correction. The same strategy was applied to analyze the differential expression of phosphoproteomic data.

Pathway enrichment analysis

Differentially expressed genes were subjected to Gene Ontology and KEGG pathway enrichment analysis in DAVID[95] with a P-value/FDR < 0.05. We used gene sets of molecular pathways from the KEGG [96] /Hallmark [97]/Reactome [98] / GO [99] databases to compute pathways.

Pathway scores and correlation analysis

Single-sample gene set enrichment analysis (ssGSEA) [100] was utilized to obtain pathway scores for each sample based on proteomic and phosphoproteomic data using the R package GSVA [101]. The KEGG, Reactome, and Hallmark gene sets from the MsigDB v7.1 were used as the background. Correlations between the pathway scores and other features were determined using Spearman’s correlation [102]. Inferred activity was performed using ssGSEA implemented in the R package GSVA with a minimum gene set size of 10.

Survival analysis

Kaplan-Meier survival curves (log-rank test) were used to determine the OS and PFS of different subtypes and patients. The coefficient value, which is equal to ln (HR), was calculated using Cox proportional hazards regression analysis. P-values less than 0.05, were considered significantly different and selected for Cox regression multivariate analysis. Prior to the log-rank test of a given protein, phosphoprotein, or phosphosite, survminer (version 0.2.4, R package) with maxstat (maximally selected rank statistics; http://r-addict.com/2016/11/21/Optimal-Cutpoint-maxstat.html) was used to determine the optimal cut-off point for the selected samples according to a previous study[103]. OS curves were then calculated (Kaplan-Meier analysis, log-rank test) based on the optimal cut-off point.

Gene Set Enrichment Analysis (GSEA)

GSEA was performed by the GSEA software (http://software.broadinstitute.org/gsea/ index.jsp). Gene sets including KEGG, GO Biological Process (BP), Reactome and HALLMARK downloaded from the Molecular Signatures Database (MSigDB v7.1, http://software.broadinstitute.org/gsea/msigdb/index.jsp) were used. FDR < 0.05 was used as a cutoff. The normalized enrichment score was used to reflect the degree of pathway overrepresentation.

Multivariate Cox proportional hazards regression analysis

Multivariate Cox proportional hazard regression analysis was used to establish the PFS-predicting model for NMIBC. The results were plotted using a forest map with the R package “ggplot2”.

Construction and validation of predictive models to distinguish between BCG-R and BCG-NR

Binomial logistic regression analysis was used to construct the distinguishing between BCG-R and BCG-NR predicting model based on the significantly differently expressed proteins in BCG-R and BCG-NR samples using in the R software v3.5.1. The backward stepwise method was utilized to feature selection. Samples were randomly divided into the training set and the testing set. Moreover, the diagnostic value of this model was verified using ROC analysis (pROC R package version 1.16.2 and Caret R package version 6.0–86). Sensitivity, specificity, accuracy, and AUC were used to determine predictive values. The predictive model was validated in validation cohort.

Phosphoproteomic data analysis

Database searching of MS phosphoproteomic data

Phospho-proteome MS raw files were searched against the human Refseq protein database (27,414 proteins, version 04/07/2013) using Proteome Discoverer (version 2.3.0.523) with a Mascot [104] (version 2.3.01) engine with a percolator [105]. Carbamidomethyl cysteine was used as a fixed modification, and oxidized methionine, protein N-term acetylation, and phospho (S/T/Y) were set as variable modifications. The FDR of peptides and proteins was set at 1%. The tolerance for spectral search a mass tolerance of 20 ppm for the precursor. The maximum number of missing cleavage site was set at 2. For phosphosite localization, ptmRS [106] was used to determine phosphosite confidence, and a phosphosite probability > 0.75 was used for further analysis.

Kinase Activity Prediction

Kinase activity scores were inferred from phosphorylation sites by employing PTM signature enrichment analysis (PTM-SEA) using the PTM signatures database (PTMsigDB) v1.9.0 (https://github.com/broadinstitute/ssGSEA2.0). Sequence windows flanking the phosphorylation site by 7 amino acids in both directions were used as unique site identifiers. Only fully localized phosphorylation sites as determined by Spectrum Mill software were taken into consideration. Phosphorylation sites on multiply phosphorylated peptides were resolved using the approach described in Krug et al.[107] resulting in a total of 29,406 phosphorylation sites that were subjected to PTM-SEA analysis using the following parameters:

gene.set.database = ‘‘ptm.sig.db.all.flanking.human.v1.9.0.gmt’’

sample.norm.type = ‘‘rank’’

weight = 0.75

statistic = ‘‘area.under.RES’’

output.score.type = ’’NES’’

nperm = 1000

global.fdr = TRUE

min.overlap = 5

correl.type = ‘‘z.score”

Immunohistochemistry (IHC)

Formalin-fixed, paraffin-embedded tissue sections of 10 µM thickness were stained in batches for detecting TLR3 in a central laboratory at the FUSCC according to standard automated protocols. Deparaffinization and rehydration were performed, followed by antigen retrieval and antibody staining. TLR3 IHC was performed using the Leica BOND-MAX auto staining system[108]. Rabbit monoclonal anti-TLR3 antibody (ABclonal, A11778) was introduced, followed by detection with a Bond Polymer Refine Detection DS9800 (Bond). Slides were imaged using an OLYMPUS BX43 microscope (OLYMPUS) and processed using a Scanscope (Leica).

Functional Experiments

Cell culture

The human bladder cancer cell lines (T24, 5637 and RT4) were purchased from American Type Culture Collection (ATCC). T24 and 5637 were cultured in Dulbecco’s modified eagle medium (DMEM, Gibco, Carlsbad, CA, United States) containing 10% fetal bovine serum (FBS, Gibco, Carlsbad, CA, United States) and 100 U/ml penicillin streptomycin. RT4 was cultured in McCoys 5A Medium with 10% fetal bovine serum (FBS, Gibco, Carlsbad, CA, United States) and 100 U/ml penicillin streptomycin. For TLR3 overexpression, the full-length TLR3 was cloned into pCDH-puro vector plasmid.

Western blot

Western blotting was performed as previously described [109]. Proteins were extracted from harvested BC cells separately and quantified by BCA assay. The primary antibodies used in this study were anti-TLR3 antibody (1:1000, CST, 6961s), anti-p65 antibody (1:1000, CST,3034), anti-p-p65 antibody (1:1000, CST,3031), HLA class I ABC Polyclonal antibody (Proteintech, 15240-1-AP) and anti-GAPDH antibody (ab181602, Abcam), β-Actin antibody (CST,4967) in our study.

Flow cytometry Analysis for HLA-I

After collecting cell precipitates and washing them once with PBS, the cells were resuspended in 500 cell staining buffer. Subsequently, they were incubated in the dark for 30 minutes with gentle agitation, then Alexa Fluor® 488 anti-human HLA-A, B, C antibody was immediately added. They were next incubated at room temperature for an additional 30 minutes. Eventually, the samples were analyzed promptly using a flow cytometer (FACSCalibur flow cytometer from BD Biosciences). Each experiment should be conducted at least three times.

Apoptosis Analysis

Apoptosis in cells was assessed using the Annexin V-PE/7-AAD Apoptosis Detection Kit (Liankebio, Hangzhou, China) and flow cytometry. Following the collection of cell precipitates, a single wash with PBS was performed. Subsequently, the cells were resuspended in 200 µl of 1x binding buffer. The resuspended cells were then incubated with Annexin V-PE and 7-AAD at room temperature for 5 minutes, followed by immediate analysis using flow cytometry.

Tissue digestion and Patient derived organoid establishment

Ten tissue specimens from BC patients undergoing transurethral cystectomy (TUR-B) at the Department of Urology, Affiliated Cancer Hospital of Fudan University, were utilized to establish Patient-Derived Organoids (PDOs). Initially, the tissue specimens were immersed in a basal medium (Advanced DMEM F12 Serum-Free medium, Gibco, 12634010). After mechanical disruption, the tumor tissues were minced into a paste, rinsed in Basal medium (800 g for 5 min), and enzymatically digested at 37°C. The digestion mixture contained 2.5 mg/ml Collagenase II (ThermoFisher 17101015) and 10 µM Y-27632-HCl Rock Inhibitor (MCE, 146986-50-7) in Basal medium for 10–15 minutes, with periodic stirring every 5 minutes. The digested tissue underwent additional washing with Basal medium (1000 rpm, 5 min). The resulting precipitate was then resuspended in 1–2 ml TrypLE Express (ThermoFisher, 12605028) and subjected to a 37°C digestion for 5 min. The digestion was terminated using Basal medium through a 70 µm cell strainer (Corning, 352350). The cell filtrate underwent centrifugation once more (1000 rpm, 5 min). The cell clusters were subsequently resuspended in an organoid matrix gel (Corning, 356231) and seeded into Ultra Low Adhesion (ULA) 24-well plates (Corning, 3743). Following a 30-minute incubation at 37°C under 5% CO2 conditions, Bladder Cancer Organoid Specific Medium (500 µl/well, Bladder Organoid Kit, K2126-CB) was added. Relative cell activity was measured using CellTiter-Glo 3D (Promega, G9682) following the manufacturer's instructions [110, 111].

ELISA

To detect IL-10, IL-4, and IFN-beta, specific ELISA kits are utilized (IFN-beta, JL19215; IL10, JL19246; IL-4, JL19870). The procedure involves coating a microplate with capture antibodies that are specific to the cytokines of interest. After coating, the plate is blocked to prevent any nonspecific binding. Next, the cell culture supernatants containing the cytokines are added to the wells, allowing the target cytokines to bind to the capture antibodies. Following a washing step to remove any unbound substances, detection antibodies conjugated with enzymes are added. These antibodies bind to the captured cytokines. After another round of washing to eliminate excess detection antibodies, a substrate solution is added to the wells. The enzyme catalyzes a reaction that produces a measurable color change. The intensity of the color is proportional to the concentration of the cytokines present in the cell culture supernatant. Finally, the absorbance is measured using a microplate reader, and the concentrations of IL-10, IL-4, and IFN-beta are determined by comparing the absorbance values to a standard curve generated with known concentrations of the respective cytokines.

In vitro BCG stimulation assay

Bladder cancer cell lines and organoids were treated with Bacillus Calmette-Guérin (BCG) at a multiplicity of infection (MOI) of 10 [49, 112]. The cells were cultured in an appropriate growth medium until they reached 70–80% confluence. A BCG solution was then prepared and added to the cells to ensure uniform exposure. Following treatment, the cells (24h) and organoids (7 days) were incubated, and after the designated period, they were harvested for analysis.

Statistical analysis

Standard statistical tests were used to analyze the clinical data, including but not limited to Student’s t test, Wilcoxon rank-sum test, Chi-square test, Fisher’s exact test, Kruskal-Wallis test, Log-rank test. For categorical variables versus categorical variables (including gene mutations, gender, age group, smoke status, nerve invasion, vascular invasion, metastasis, hyperglycemia, hypertension, TNM stage, and histological type), Fisher’s exact test was used in a 2 × 2 table, otherwise Chi-square test was used. The Wilcoxon rank-sum test was used to examine whether genes were differentially expressed between different proteomic subtypes, and patients with different mutation statuses and CNA of statuses. The Kruskal-Wallis test was used to test whether genes were differentially expressed among the different tissues type or other subgroups. To account for multiple-testing, the P values were adjusted using the Benjamini-Hochberg FDR correction. Kaplan–Meier plots (Log-rank test) were used to describe OS. Variables associated with OS and PFS were identified using univariate Cox proportional hazards regression models. Significant factors in univariate analysis were further subjected to a multivariate Cox regression analysis in a forward LR manner. All the analyses of clinical data were performed in R (version 3.5.1). For functional experiments, three biological repeats were performed independently, and results were expressed as mean ± standard error of the mean (SEM). The statistical significance of differences was determined by two-way ANOVA. Statistical analysis was performed using GraphPad Prism (version 5.01). The p values less than 0.05, 0.01, 0.001, 0.0001 were marked with *, **, ***, ****, respectively. All the statistical analyses had been checked by two statisticians.

FFPE

Formalin-fixed paraffin-embedded

FOT

Fraction of total

FDR

False discovery rate

CNA

Copy number alteration

FPKM

Fragments Per Kilobase of transcript per Million mapped reads

HPLC

High-performance liquid chromatography

iBAQ

Intensity-based absolute quantification

IHC

Immunohistochemistry

Mass spectrometry

PCA

Principal-component analysis

ssGSEA

Single sample gene set enrichment analysis

Urothelial cancer

Ethics approval and consent to participate

The present study was carried out comply with the ethical standards of Helsinki Declaration II and

approved by the Institution Review Board of FUSCC (050432-4-1911D).

Availability of data and material

Proteome and phosphoprotome raw datasets have been deposited to the ProteomeXchange Consortium (Project ID: IPX0006877000) via the iProX partner repository (https://www.iprox.cn). WES data were deposited in GSA-human (Genome Sequence Archive for human) in NGDC (the National Genomics Data Center) (https://ngdc.cncb.ac.cn/). The raw sequencing data are available under controlled access due to data privacy laws related to patient consent for data sharing and the data should be used for research purposes only. Access can be obtained by approval via their respective DAC (Data Access Committees) in the GSA-human database. According to the guidelines of GSA-human, all non-profit researchers are allowed access to the data and the Principle Investigator of any research group is allowed to apply for Controlled access of the data. The user can register and login to the GSA database website (https://ngdc.cncb.ac.cn/gsa-human/) and follow the guidance of “Request Data” to request the data step by step (https://ngdc.cncb.ac.cn/gsa-human/document/GSA-Human_Request_Guide_for_Users_us.pdf). The approximate response time for accession requests is about 2 weeks. The access authority can be obtained for Research Use Only. The user can also contact the corresponding author directly. Once access has been granted, the data will be available to download for 3 months. The remaining data are available within the Article, Supplementary Information, or Source Data file. Source data are provided with this paper.

Declaration of interests

The authors declare that they have no competing interests.

Acknowledgements

This work is supported by the National Key Research and Development Program of China (2022YFA1303200 [C.D.], 2022YFA1303201 [C.D.]), the National Natural Science Foundation of China: 32330062 [C.D.], 31972933 [C.D.], 82172817 [Y.Y.Q.], 82172741 [Y.Y.Q.], the National Ten Thousand Plan Young Top Talents [Y.Y.Q.]，Program of Shanghai Academic/Technology Research Leader (22XD1420100) [C.D.], the Major Project of Special Development Funds of Zhangjiang National Independent Innovation Demonstration Zone (ZJ2019‐ZD‐004 [C.D.]), the Science and Technology Commission of Shanghai Municipality (2017SHZDZX01 [C.D.]), the Young Scientists Fund of the National Natural Science Foundation of China (32201215 [J.W.F.] and 32201212 [Y.Z.W.], Shanghai Rising-Star Program (23QA1408900) [Y.Y.Q.], Shanghai "Science and Technology Innovation Action Plan" medical innovation research Project (22Y11905100) [Y.Y.Q.], Shanghai Municipal Health Bureau Project (No. 2020CXJQ03) [Y.Y.Q.]. Shanghai Sailing Program (22YF1403100 [J.W.F.]), the China Postdoctoral Science Foundation:2023TQ0084 [N.X.], and the Fudan original research personalized support project.

Authors’ information

Department of Urology, Fudan University Shanghai Cancer Center, State Key Laboratory of Genetic Engineering, Collaborative Innovation Center for Genetics and Development, School of Life Sciences, Institutes of Biomedical Sciences, and Human Phenome Institute, Fudan University, Shanghai 20433, China; Department of Oncology, Shanghai Medical College, Fudan University, Shanghai Genitourinary Cancer Institute, Shanghai 200032, China; Tissue Bank & Department of Pathology, Fudan University Shanghai Cancer Center, Shanghai 200032, China; Institute for Development and Regenerative Cardiovascular Medicine, MOE-Shanghai Key Laboratory of Children’s Environmental Health, Xinhua Hospital, Shanghai Jiao Tong University School of Medicine, Shanghai 200092, China.

Corresponding authors

Correspondence to Yuanyuan Qu, Jianyuan Zhao, Dingwei Ye or Chen Ding.

Contributions

Conceptualization: D.Y., C.D., J.G., J.Z., and Y.Q., Performed Experiment and Data Collection: Y.Q., N.X., Z.Y., T. C., and Y.W., Data Curation: J.F., Y.W., N.X., Z.Y., Y.Q., and C.D., Data Analysis: Z.Y., N.X., T.C., L.Z., and Y.W., Visualization: Y.Q., N.X., Z.Y., T.C., Y.W., Z.Q., J.Z., Y.Q., X.X., S.T and Q.Z., Patient Sample Management and QC: L.Z., Y.Z., C.G., Y.S., Y.G., S.N., and F.D., Supervision: H.Z., C.D., J.Z., and D.Y., Writing: Z.Y., N.X., L.Z., Y.Q., and C.D.

Teoh JY, Huang J, Ko WY, Lok V, Choi P, Ng CF, Sengupta S, Mostafid H, Kamat AM, Black PC, et al. Global Trends of Bladder Cancer Incidence and Mortality, and Their Associations with Tobacco Use and Gross Domestic Product Per Capita. Eur Urol. 2020;78(6):893–906.
Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, Bray F. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin. 2021;71(3):209–49.
Antoni S, Ferlay J, Soerjomataram I, Znaor A, Jemal A, Bray F. Bladder Cancer Incidence and Mortality: A Global Overview and Recent Trends. Eur Urol. 2017;71(1):96–108.
DeGeorge KC, Holt HR, Hodges SC. Bladder Cancer: Diagnosis and Treatment. Am Fam Physician. 2017;96(8):507–14.
Ansari Djafari A, Javanmard B, Razzaghi M, Hojjati SA, Razzaghi Z, Faraji S, Rahavian A, Garoosi M. Intravesical Gemcitabine versus Intravesical Bacillus Calmette-Guerin for the Treatment of Intermediate-Risk Non-Muscle Invasive Bladder Cancer: A Randomized Controlled Trial. Urol J. 2023;20(2):123–8.
van den Bosch S, Alfred Witjes J. Long-term cancer-specific survival in patients with high-risk, non-muscle-invasive bladder cancer and tumour progression: a systematic review. Eur Urol. 2011;60(3):493–500.
Ritch CR, Velasquez MC, Kwon D, Becerra MF, Soodana-Prakash N, Atluri VS, Almengo K, Alameddine M, Kineish O, Kava BR, et al. Use and Validation of the AUA/SUO Risk Grouping for Nonmuscle Invasive Bladder Cancer in a Contemporary Cohort. J Urol. 2020;203(3):505–11.
Sylvester RJ, van der Lamm MA. Intravesical bacillus Calmette-Guerin reduces the risk of progression in patients with superficial bladder cancer: a meta-analysis of the published results of randomized clinical trials. J Urol. 2002;168(5):1964–70.
Sylvester RJ, van der Meijden AP, Witjes JA, Kurth K. Bacillus calmette-guerin versus chemotherapy for the intravesical treatment of patients with carcinoma in situ of the bladder: a meta-analysis of the published results of randomized clinical trials. J Urol. 2005;174(1):86–91. discussion 91 – 82.
Pettenati C, Ingersoll MA. Mechanisms of BCG immunotherapy and its outlook for bladder cancer. Nat Rev Urol. 2018;15(10):615–25.
Bajic P, Wolfe AJ, Gupta GN. The Urinary Microbiome: Implications in Bladder Cancer Pathogenesis and Therapeutics. Urology. 2019;126:10–5.
Vandeveer AJ, Fallon JK, Tighe R, Sabzevari H, Schlom J, Greiner JW. Systemic Immunotherapy of Non-Muscle Invasive Mouse Bladder Cancer with Avelumab, an Anti-PD-L1 Immune Checkpoint Inhibitor. Cancer Immunol Res. 2016;4(5):452–62.
Rentsch CA, Derre L, Dugas SG, Wetterauer C, Federer-Gsponer JR, Thalmann GN, Ingersoll MA. Building on a Solid Foundation: Enhancing Bacillus Calmette-Guerin Therapy. Eur Urol Focus. 2018;4(4):485–93.
Witjes JA. Management of BCG failures in superficial bladder cancer: a review. Eur Urol. 2006;49(5):790–7.
Tse J, Singla N, Ghandour R, Lotan Y, Margulis V. Current advances in BCG-unresponsive non-muscle invasive bladder cancer. Expert Opin Investig Drugs. 2019;28(9):757–70.
Fernandez-Gomez J, Madero R, Solsona E, Unda M, Martinez-Pineiro L, Gonzalez M, Portillo J, Ojea A, Pertusa C, Rodriguez-Molina J, et al. Predicting nonmuscle invasive bladder cancer recurrence and progression in patients treated with bacillus Calmette-Guerin: the CUETO scoring model. J Urol. 2009;182(5):2195–203.
Lindskrog SV, Prip F, Lamy P, Taber A, Groeneveld CS, Birkenkamp-Demtroder K, Jensen JB, Strandgaard T, Nordentoft I, Christensen E, et al. An integrated multi-omics analysis identifies prognostic molecular subtypes of non-muscle-invasive bladder cancer. Nat Commun. 2021;12(1):2301.
Robertson AG, Groeneveld CS, Jordan B, Lin X, McLaughlin KA, Das A, Fall LA, Fantini D, Taxter TJ, Mogil LS, et al. Identification of Differential Tumor Subtypes of T1 Bladder Cancer. Eur Urol. 2020;78(4):533–7.
Babjuk M, Burger M, Capoun O, Cohen D, Comperat EM, Dominguez Escrig JL, Gontero P, Liedberg F, Masson-Lecomte A, Mostafid AH, et al. European Association of Urology Guidelines on Non-muscle-invasive Bladder Cancer (Ta, T1, and Carcinoma in Situ). Eur Urol. 2022;81(1):75–94.
West AP, Koblansky AA, Ghosh S. Recognition and signaling by toll-like receptors. Annu Rev Cell Dev Biol. 2006;22:409–37.
Pietzak EJ, Bagrodia A, Cha EK, Drill EN, Iyer G, Isharwal S, Ostrovnaya I, Baez P, Li Q, Berger MF, et al. Next-generation Sequencing of Nonmuscle Invasive Bladder Cancer Reveals Potential Biomarkers and Rational Therapeutic Targets. Eur Urol. 2017;72(6):952–9.
Levatic J, Salvadores M, Fuster-Tormo F, Supek F. Mutational signatures are markers of drug sensitivity of cancer cells. Nat Commun. 2022;13(1):2926.
Samstein RM, Lee CH, Shoushtari AN, Hellmann MD, Shen R, Janjigian YY, Barron DA, Zehir A, Jordan EJ, Omuro A, et al. Tumor mutational load predicts survival after immunotherapy across multiple cancer types. Nat Genet. 2019;51(2):202–6.
Xu N, Yao Z, Shang G, Ye D, Wang H, Zhang H, Qu Y, Xu F, Wang Y, Qin Z, et al. Integrated proteogenomic characterization of urothelial carcinoma of the bladder. J Hematol Oncol. 2022;15(1):76.
Haura EB, Turkson J, Jove R. Mechanisms of disease: Insights into the emerging role of signal transducers and activators of transcription in cancer. Nat Clin Pract Oncol. 2005;2(6):315–24.
Verhoeven Y, Tilborghs S, Jacobs J, De Waele J, Quatannens D, Deben C, Prenen H, Pauwels P, Trinh XB, Wouters A, et al. The potential and controversy of targeting STAT family members in cancer. Semin Cancer Biol. 2020;60:41–56.
Ferrington DA, Gregerson DS. Immunoproteasomes: structure, function, and antigen presentation. Prog Mol Biol Transl Sci. 2012;109:75–112.
Tanaka N, Taniguchi T. The interferon regulatory factors and oncogenesis. Semin Cancer Biol. 2000;10(2):73–81.
Brucet M, Marques L, Sebastian C, Lloberas J, Celada A. Regulation of murine Tap1 and Lmp2 genes in macrophages by interferon gamma is mediated by STAT1 and IRF-1. Genes Immun. 2004;5(1):26–35.
Hayashi T, Horiuchi A, Sano K, Hiraoka N, Kasai M, Ichimura T, Sudo T, Tagawa Y, Nishimura R, Ishiko O, et al. Potential role of LMP2 as tumor-suppressor defines new targets for uterine leiomyosarcoma therapy. Sci Rep. 2011;1:180.
Hurst CD, Platt FM, Taylor CF, Knowles MA. Novel tumor subgroups of urothelial carcinoma of the bladder defined by integrated genomic analysis. Clin Cancer Res. 2012;18(21):5865–77.
Iwasaki A, Medzhitov R. Toll-like receptor control of the adaptive immune responses. Nat Immunol. 2004;5(10):987–95.
Ayari C, Bergeron A, LaRue H, Menard C, Fradet Y. Toll-like receptors in normal and malignant human bladders. J Urol. 2011;185(5):1915–21.
Meylan E, Burns K, Hofmann K, Blancheteau V, Martinon F, Kelliher M, Tschopp J. RIP1 is an essential mediator of Toll-like receptor 3-induced NF-kappa B activation. Nat Immunol. 2004;5(5):503–7.
Zheng X, Li S, Yang H. Roles of Toll-Like Receptor 3 in Human Tumors. Front Immunol. 2021;12:667454.
Oeckinghaus A, Hayden MS, Ghosh S. Crosstalk in NF-kappaB signaling pathways. Nat Immunol. 2011;12(8):695–708.
Hayden MS, Ghosh S. Signaling to NF-kappaB. Genes Dev. 2004;18(18):2195–224.
Israel A. The IKK complex, a central regulator of NF-kappaB activation. Cold Spring Harb Perspect Biol. 2010;2(3):a000158.
Robertson AG, Groeneveld CS, Jordan B, Lin X, McLaughlin KA, Das A, Fall LA, Fantini D, Taxter TJ, Mogil LS, et al. Identification of Differential Tumor Subtypes of T1 Bladder Cancer. Eur Urol. 2020;78(4):533–7.
Lindskrog SV, Prip F, Lamy P, Taber A, Groeneveld CS, Birkenkamp-Demtröder K, Jensen JB, Strandgaard T, Nordentoft I, Christensen E, et al. An integrated multi-omics analysis identifies prognostic molecular subtypes of non-muscle-invasive bladder cancer. Nat Commun. 2021;12(1):2301.
Damrauer JS, Hoadley KA, Chism DD, Fan C, Tiganelli CJ, Wobker SE, Yeh JJ, Milowsky MI, Iyer G, Parker JS, et al. Intrinsic subtypes of high-grade bladder cancer reflect the hallmarks of breast cancer biology. Proc Natl Acad Sci U S A. 2014;111(8):3110–5.
Kamoun A, de Reyniès A, Allory Y, Sjödahl G, Robertson AG, Seiler R, Hoadley KA, Groeneveld CS, Al-Ahmadie H, Choi W, et al. A Consensus Molecular Classification of Muscle-invasive Bladder Cancer. Eur Urol. 2020;77(4):420–33.
Xing X, Cai L, Ouyang J, Wang F, Li Z, Liu M, Wang Y, Zhou Y, Hu E, Huang C, et al. Proteomics-driven noninvasive screening of circulating serum protein panels for the early diagnosis of hepatocellular carcinoma. Nat Commun. 2023;14(1):8392.
Hoshino A, Kim HS, Bojmar L, Gyan KE, Cioffi M, Hernandez J, Zambirinis CP, Rodrigues G, Molina H, Heissel S, et al. Extracellular Vesicle and Particle Biomarkers Define Multiple Human Cancers. Cell. 2020;182(4):1044–1061e1018.
Liu Y, Gu Y, Han Y, Zhang Q, Jiang Z, Zhang X, Huang B, Xu X, Zheng J, Cao X. Tumor Exosomal RNAs Promote Lung Pre-metastatic Niche Formation by Activating Alveolar Epithelial TLR3 to Recruit Neutrophils. Cancer Cell. 2016;30(2):243–56.
Wu W, Chen L, Jia G, Tang Q, Han B, Xia S, Jiang Q, Liu H. Inhibition of FGFR3 upregulates MHC-I and PD-L1 via TLR3/NF-kB pathway in muscle-invasive bladder cancer. Cancer Med. 2023;12(14):15676–90.
El-Shitany NA, Eid BG. Icariin modulates carrageenan-induced acute inflammation through HO-1/Nrf2 and NF-kB signaling pathways. Biomed Pharmacother. 2019;120:109567.
Au L, Hatipoglu E, Robert de Massy M, Litchfield K, Beattie G, Rowan A, Schnidrig D, Thompson R, Byrne F, Horswell S, et al. Determinants of anti-PD-1 response and resistance in clear cell renal cell carcinoma. Cancer Cell. 2021;39(11):1497–1518e1411.
Rouanne M, Adam J, Radulescu C, Letourneur D, Bredel D, Mouraud S, Goubet AG, Leduc M, Chen N, Tan TZ et al. BCG therapy downregulates HLA-I on malignant cells to subvert antitumor immune responses in bladder cancer. J Clin Invest 2022, 132(12).
Ouyang W, Rutz S, Crellin NK, Valdez PA, Hymowitz SG. Regulation and functions of the IL-10 family of cytokines in inflammation and disease. Annu Rev Immunol. 2011;29:71–109.
Chen Y, Lin J, Zhao Y, Ma X, Yi H. Toll-like receptor 3 (TLR3) regulation mechanisms and roles in antiviral innate immune responses. J Zhejiang Univ Sci B. 2021;22(8):609–32.
Bernardo AR, Cosgaya JM, Aranda A, Jimenez-Lara AM. Synergy between RA and TLR3 promotes type I IFN-dependent apoptosis through upregulation of TRAIL pathway in breast cancer cells. Cell Death Dis. 2013;4(1):e479.
Muthuswamy R, Wang L, Pitteroff J, Gingrich JR, Kalinski P. Combination of IFNalpha and poly-I:C reprograms bladder cancer microenvironment for enhanced CTL attraction. J Immunother Cancer. 2015;3:6.
Rizvi NA, Hellmann MD, Snyder A, Kvistborg P, Makarov V, Havel JJ, Lee W, Yuan J, Wong P, Ho TS, et al. Cancer immunology. Mutational landscape determines sensitivity to PD-1 blockade in non-small cell lung cancer. Science. 2015;348(6230):124–8.
Carbone DP, Reck M, Paz-Ares L, Creelan B, Horn L, Steins M, Felip E, van den Heuvel MM, Ciuleanu TE, Badin F, et al. First-Line Nivolumab in Stage IV or Recurrent Non-Small-Cell Lung Cancer. N Engl J Med. 2017;376(25):2415–26.
Yarchoan M, Hopkins A, Jaffee EM. Tumor Mutational Burden and Response Rate to PD-1 Inhibition. N Engl J Med. 2017;377(25):2500–1.
Hellmann MD, Callahan MK, Awad MM, Calvo E, Ascierto PA, Atmaca A, Rizvi NA, Hirsch FR, Selvaggi G, Szustakowski JD, et al. Tumor Mutational Burden and Efficacy of Nivolumab Monotherapy and in Combination with Ipilimumab in Small-Cell Lung Cancer. Cancer Cell. 2018;33(5):853–861e854.
Hellmann MD, Ciuleanu TE, Pluzanski A, Lee JS, Otterson GA, Audigier-Valette C, Minenza E, Linardou H, Burgers S, Salman P, et al. Nivolumab plus Ipilimumab in Lung Cancer with a High Tumor Mutational Burden. N Engl J Med. 2018;378(22):2093–104.
Gobin SJ, Peijnenburg A, Keijsers V, van den Elsen PJ. Site alpha is crucial for two routes of IFN gamma-induced MHC class I transactivation: the ISRE-mediated route and a novel pathway involving CIITA. Immunity. 1997;6(5):601–11.
Solheim JC. Class I MHC molecules: assembly and antigen presentation. Immunol Rev. 1999;172:11–9.
Chatterjee-Kishore M, Wright KL, Ting JP, Stark GR. How Stat1 mediates constitutive gene expression: a complex of unphosphorylated Stat1 and IRF1 supports transcription of the LMP2 gene. EMBO J. 2000;19(15):4111–22.
Osterloh P, Linkemann K, Tenzer S, Rammensee HG, Radsak MP, Busch DH, Schild H. Proteasomes shape the repertoire of T cells participating in antigen-specific immune responses. Proc Natl Acad Sci U S A. 2006;103(13):5042–7.
Groettrup M, Kirk CJ, Basler M. Proteasomes in immune cells: more than peptide producers? Nat Rev Immunol. 2010;10(1):73–8.
Basler M, Beck U, Kirk CJ, Groettrup M. The antiviral immune response in mice devoid of immunoproteasome activity. J Immunol. 2011;187(11):5548–57.
Ayari C, Besancon M, Bergeron A, LaRue H, Bussieres V, Fradet Y. Poly(I:C) potentiates Bacillus Calmette-Guerin immunotherapy for bladder cancer. Cancer Immunol Immunother. 2016;65(2):223–34.
Cebrian I, Croce C, Guerrero NA, Blanchard N, Mayorga LS. Rab22a controls MHC-I intracellular trafficking and antigen cross-presentation by dendritic cells. EMBO Rep. 2016;17(12):1753–65.
Wu SJ, Niknafs YS, Kim SH, Oravecz-Wilson K, Zajac C, Toubai T, Sun Y, Prasad J, Peltier D, Fujiwara H, et al. A Critical Analysis of the Role of SNARE Protein SEC22B in Antigen Cross-Presentation. Cell Rep. 2017;19(13):2645–56.
Elhamamsy AR, Metge BJ, Alsheikh HA, Shevde LA, Samant RS. Ribosome Biogenesis: A Central Player in Cancer Metastasis and Therapeutic Resistance. Cancer Res. 2022;82(13):2344–53.
Jiao L, Liu Y, Yu XY, Pan X, Zhang Y, Tu J, Song YH, Li Y. Ribosome biogenesis in disease: new players and therapeutic targets. Signal Transduct Target Ther. 2023;8(1):15.
Chamie K, Litwin MS, Bassett JC, Daskivich TJ, Lai J, Hanley JM, Konety BR, Saigal CS. Urologic Diseases in America P: Recurrence of high-risk bladder cancer: a population-based analysis. Cancer. 2013;119(17):3219–27.
de Jong FC, Laajala TD, Hoedemaeker RF, Jordan KR, van der Made ACJ, Boeve ER, van der Schoot DKE, Nieuwkamer B, Janssen EAM, Mahmoudi T, et al. Non-muscle-invasive bladder cancer molecular subtypes predict differential response to intravesical Bacillus Calmette-Guerin. Sci Transl Med. 2023;15(697):eabn4118.
Lerner SP, Dinney C, Kamat A, Bivalacqua TJ, Nielsen M, O'Donnell M, Schoenberg MP, Steinberg G. Clarification of Bladder Cancer Disease States Following Treatment of Patients with Intravesical BCG. Bladder Cancer. 2015;1(1):29–30.
Kamat AM, Sylvester RJ, Bohle A, Palou J, Lamm DL, Brausi M, Soloway M, Persad R, Buckley R, Colombel M, et al. Definitions, End Points, and Clinical Trial Designs for Non-Muscle-Invasive Bladder Cancer: Recommendations From the International Bladder Cancer Group. J Clin Oncol. 2016;34(16):1935–44.
Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009;25(14):1754–60.
Etherington GJ, Ramirez-Gonzalez RH, MacLean D. bio-samtools 2: a package for analysis and visualization of sequence and alignment data with SAMtools in Ruby. Bioinformatics. 2015;31(15):2565–7.
Koboldt DC, Zhang Q, Larson DE, Shen D, McLellan MD, Lin L, Miller CA, Mardis ER, Ding L, Wilson RK. VarScan 2: somatic mutation and copy number alteration discovery in cancer by exome sequencing. Genome Res. 2012;22(3):568–76.
Cibulskis K, Lawrence MS, Carter SL, Sivachenko A, Jaffe D, Sougnez C, Gabriel S, Meyerson M, Lander ES, Getz G. Sensitive detection of somatic point mutations in impure and heterogeneous cancer samples. Nat Biotechnol. 2013;31(3):213–9.
Sherry ST, Ward MH, Kholodov M, Baker J, Phan L, Smigielski EM, Sirotkin K. dbSNP: the NCBI database of genetic variation. Nucleic Acids Res. 2001;29(1):308–11.
Abecasis GR, Auton A, Brooks LD, DePristo MA, Durbin RM, Handsaker RE, Kang HM, Marth GT, McVean GA. An integrated map of genetic variation from 1,092 human genomes. Nature. 2012;491(7422):56–65.
Adzhubei I, Jordan DM, Sunyaev SR. Predicting functional effect of human missense mutations using PolyPhen-2. Curr Protoc Hum Genet 2013, Chap. 7:Unit7.20.
Kircher M, Witten DM, Jain P, O'Roak BJ, Cooper GM, Shendure J. A general framework for estimating the relative pathogenicity of human genetic variants. Nat Genet. 2014;46(3):310–5.
Mermel CH, Schumacher SE, Hill B, Meyerson ML, Beroukhim R, Getz G. GISTIC2.0 facilitates sensitive and confident localization of the targets of focal somatic copy-number alteration in human cancers. Genome Biol. 2011;12(4):R41.
Vasaikar S, Huang C, Wang X, Petyuk VA, Savage SR, Wen B, Dou Y, Zhang Y, Shi Z, Arshad OA, et al. Proteogenomic Analysis of Human Colon Cancer Reveals New Therapeutic Opportunities. Cell. 2019;177(4):1035–1049e1019.
Lawrence MS, Stojanov P, Polak P, Kryukov GV, Cibulskis K, Sivachenko A, Carter SL, Stewart C, Mermel CH, Roberts SA, et al. Mutational heterogeneity in cancer and the search for new cancer-associated genes. Nature. 2013;499(7457):214–8.
Wang S, Tao Z, Wu T, Liu XS. Sigflow: an automated and comprehensive pipeline for cancer genome mutational signature analysis. Bioinformatics. 2021;37(11):1590–2.
Rosenthal R, McGranahan N, Herrero J, Taylor BS, Swanton C. DeconstructSigs: delineating mutational processes in single tumors distinguishes DNA repair deficiencies and patterns of carcinoma evolution. Genome Biol. 2016;17(1):1–11.
Wang K, Li M, Hakonarson H. ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data. Nucleic Acids Res. 2010;38(16):e164.
Chalmers ZR, Connelly CF, Fabrizio D, Gay L, Ali SM, Ennis R, Schrock A, Campbell B, Shlien A, Chmielecki J, et al. Analysis of 100,000 human cancer genomes reveals the landscape of tumor mutational burden. Genome Med. 2017;9(1):34.
Wisniewski JR, Zougman A, Nagaraj N, Mann M. Universal sample preparation method for proteome analysis. Nat Methods. 2009;6(5):359–62.
Gao Q, Zhu H, Dong L, Shi W, Chen R, Song Z, Huang C, Li J, Dong X, Zhou Y, et al. Integrated Proteogenomic Characterization of HBV-Related Hepatocellular Carcinoma. Cell. 2019;179(5):1240.
Feng J, Ding C, Qiu N, Ni X, Zhan D, Liu W, Xia X, Li P, Lu B, Zhao Q. Firmiana: towards a one-stop proteomic cloud platform for data processing and analysis. Nat Biotechnol. 2017;35(5):409–12.
Tyanova S, Temu T, Cox J. The MaxQuant computational platform for mass spectrometry-based shotgun proteomics. Nat Protoc. 2016;11(12):2301–19.
Ge S, Xia X, Ding C, Zhen B, Zhou Q, Feng J, Yuan J, Chen R, Li Y, Ge Z, et al. A proteomic landscape of diffuse-type gastric cancer. Nat Commun. 2018;9(1):1012.
Jiang Y, Sun A, Zhao Y, Ying W, Sun H, Yang X, Xing B, Sun W, Ren L, Hu B, et al. Proteomics identifies new therapeutic targets of early-stage hepatocellular carcinoma. Nature. 2019;567(7747):257–61.
Huang da W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4(1):44–57.
Ogata H, Goto S, Sato K, Fujibuchi W, Bono H, Kanehisa M. KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res. 1999;27(1):29–34.
Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, Tamayo P. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst. 2015;1(6):417–25.
Croft D, Mundo AF, Haw R, Milacic M, Weiser J, Wu G, Caudy M, Garapati P, Gillespie M, Kamdar MR, et al. The Reactome pathway knowledgebase. Nucleic Acids Res. 2014;42(Database issue):D472–477.
Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT, et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet. 2000;25(1):25–9.
Barbie DA, Tamayo P, Boehm JS, Kim SY, Moody SE, Dunn IF, Schinzel AC, Sandy P, Meylan E, Scholl C, et al. Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature. 2009;462(7269):108–12.
Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013;14:7.
Liu L, Zhu S. Computational Methods for Prediction of Human Protein-Phenotype Associations: A Review. Phenomics. 2021;1(4):171–85.
Seckinger A, Meissner T, Moreaux J, Depeweg D, Hillengass J, Hose K, Rème T, Rösen-Wolff A, Jauch A, Schnettler R, et al. Clinical and prognostic role of annexin A2 in multiple myeloma. Blood. 2012;120(5):1087–94.
Perkins DN, Pappin DJ, Creasy DM, Cottrell JS. Probability-based protein identification by searching sequence databases using mass spectrometry data. Electrophoresis. 1999;20(18):3551–67.
Käll L, Canterbury JD, Weston J, Noble WS, MacCoss MJ. Semi-supervised learning for peptide identification from shotgun proteomics datasets. Nat Methods. 2007;4(11):923–5.
Taus T, Köcher T, Pichler P, Paschke C, Schmidt A, Henrich C, Mechtler K. Universal and confident phosphorylation site localization using phosphoRS. J Proteome Res. 2011;10(12):5354–62.
Krug K, Mertins P, Zhang B, Hornbeck P, Raju R, Ahmad R, Szucs M, Mundt F, Forestier D, Jane-Valbuena J, et al. A Curated Resource for Phosphosite-specific Signature Analysis. Mol Cell Proteom. 2019;18(3):576–93.
Rabelo-Fernández RJ, Santiago-Sánchez GS, Sharma RK, Roche-Lima A, Carrion KC, Rivera RAN, Quiñones-Díaz BI, Rajasekaran S, Siddiqui J, Miles W et al. Reduced RBPMS Levels Promote Cell Proliferation and Decrease Cisplatin Sensitivity in Ovarian Cancer Cells. Int J Mol Sci 2022, 23(1).
Pillai-Kastoori L, Schutz-Geschwender AR, Harford JA. A systematic approach to quantitative Western blot analysis. Anal Biochem. 2020;593:113608.
Cattaneo CM, Dijkstra KK, Fanchi LF, Kelderman S, Kaing S, van Rooij N, van den Brink S, Schumacher TN, Voest EE. Tumor organoid-T-cell coculture systems. Nat Protoc. 2020;15(1):15–39.
Minoli M, Cantore T, Hanhart D, Kiener M, Fedrizzi T, La Manna F, Karkampouna S, Chouvardas P, Genitsch V, Rodriguez-Calero A, et al. Bladder cancer organoids as a functional system to model different disease stages and therapy response. Nat Commun. 2023;14(1):2214.
Liu K, Sun E, Lei M, Li L, Gao J, Nian X, Wang L. BCG-induced formation of neutrophil extracellular traps play an important role in bladder cancer treatment. Clin Immunol. 2019;201:4–14.
Sperber H, Mathieu J, Wang Y, Ferreccio A, Hesson J, Xu Z, Fischer KA, Devi A, Detraux D, Gu H, et al. The metabolome regulates the epigenetic landscape during naive-to-primed human embryonic stem cell transition. Nat Cell Biol. 2015;17(12):1523–35.

Table 1 is available in the Supplementary Files section.

No competing interests reported.

Table1.xlsx
TableS10119.xlsx
Supplemental table 1. Clinicopathoglogic information and multi-omics data matrix of NMIC cohort. Related to Figure 1 Related to Figure 1, Supplemental Figure 1. A. Clinical information of this cohort. B. Information of mutations identified by WES in this cohort. C. Sequencing depth of tumors with BCG-R and BCG-NR. D. Genome signature analysis of NMIBC cohort. E. Genes identified at protein level (unique peptide ≥2). F. Phosphosites identified in NMIBC samples. G. HEK293T used for QC.
TableS20116.xlsx
Supplemental table 2. The information of TMB relevant to drug sensitivity in BCG cohort. Related to Figure 2, Supplemental Figure 2. A. The TMB information in this cohort. B. The spearman correlation between TMB and enriched pathways. C. The spearman correlation between TMB and STAT1 activity. D. List of target genes of STAT1.
TableS30104.xlsx
Supplemental table 3. CNV information in BCG cohort. Related to Figure 3, Supplemental Figure 3. A. This sheet contains somatic copy number gains and the relevant genes. B. This sheet includes somatic copy number losses and the relevant genes. C. List of arm level somatic copy-number alternations events. D. List of KEGG pathways correlated with the abundance of TLR3 protein.
TableS40104.xlsx
Supplemental table 4. Multi-omics characteristics of proteomic subtypes in BCG cohort. Related to Figure 4, Supplemental Figure 4. A. Information of clusters. B. Lists of up-regulated proteins in each proteomic subtype. C. Lists of MHC-I and MHC-II proteins in the proteomic subtypes.
TableS50104.xlsx
Supplemental table 5. Differentially expressed proteins between BCG Non-Responder tumors and BCG Responder tumors. Related to Figure 5, Supplemental Figure 5. A. List of differentially proteins (Wilcoxon rank-sum test, FDR<0.05) in BCG NR tumors and BCG R tumors. B. List of candidate proteins of classifier model. C. A list of targeted peptides that unique to the proteins of classifier model.
TableS60116.xlsx
Supplemental table 6. Progression Clock of patients with NMIBC received BCG therapy. Related to Figure 6. A. The list of genes from each module. B. List of candidate proteins of the predicted model. C. A list of targeted peptides that unique to the proteins of the predicted model.
BCG.pdf
Supplementaryfigures.docx

Download PDF

Version 1

posted

You are reading this latest preprint version

Proteogenomic characterization of the non-muscle-invasive bladder cancer response to BCG reveals potential therapeutic strategies

Status:

Version 1

Abstract

Background

Methods

Results

Conclusions

Figures

Background

Results

High mutational burden (TMB) is relevant to drug sensitivity

The expression level of TLR3 in tumors is associated with prognosis and BCG treatment response

Progression Clock of Patients with NMIBC received BCG therapy

TLR3-Mediated Regulation of Cytokine Secretion in Bladder Cancer Cell Lines

Enhanced Sensitivity of TLR-High Expressing Organoids to BCG

Discussion

Materials and Methods

Whole-exome Sequencing

WES data analysis

Proteomic and phosphoproteomic analysis

MS database searching

Quantification and Statistical Analysis

Phosphoproteomic data analysis

Functional Experiments

Statistical analysis

Abbreviations

Declarations

References

Table 1

Additional Declarations

Supplementary Files

Status:

Version 1