Construction and experimental validation of a B cell-related gene signature to predict the prognosis and immunotherapeutic sensitivity in bladder cancer

Data collection and processing

The transcriptome sequencing of 19 porocarcinoma and 411 BLCA samples in “count” and “FPKM” format, along with the corresponding clinical and follow-up information, was obtained from TCGA (https://portal.gdc.cancer.gov/), as the training dataset. The gene expression matrices with the “count” format and the “FPKM” format were utilized for gene expression difference detection and prognostic model construction, respectively. We queried GEO (https://www.ncbi.nlm.nih.gov/geo/) to download the external validation datasets using the keyword “bladder cancer.” The datasets would be included if they met the following criteria: (1) The transcriptome sequencing data of the patients with BLCA was publicly available. (2) The survival statuses and follow-up duration or their clinical response to immunotherapy can be downloaded in the GEO database in the Supplementary Materials of the original article. (3) The expressions of the genes in the established model must be included in the dataset. Two researchers (Zhou R and Zhou J) independently extracted the dataset and reached a consensus for all items, and any conflicts that arose were resolved with an expert (Tan W) invited to join the discussion.

After the careful and manual review of retrieval results in the GEO, four datasets, including GSE13507 [11, 12], GSE31684 [13, 14], GSE32894 [15], and GSE111636, were extracted. The detailed information on these GEO datasets is displayed in Supplementary Table 1. Given that the GSE111636 dataset, which included 6 responders and 5 non-responders to pembrolizumab, has not been publicly published, the IMvigor210 cohort, containing 195 BLCA samples, was used in this study, whose transcription data and clinical information can be obtained from the IMvigor210CoreBiologies package in R software (version 4.0.3) [16]. In these cohorts, the cases with less than 30 days’ follow-up and the gene with average expression values < 0.5 were excluded.

Clinical sample collection

The tumor and the adjacent normal tissues from 19 BLCA patients undergoing partial/radical cystectomy without prior radiotherapy and chemotherapy were collected from the Nanfang Hospital of Southern Medical University (Guangzhou, China) between 2020 and 2022, which was named “NH cohort.” The BLCA samples were isolated from the center of the tumor and immediately stored in liquid nitrogen. The formalin-fixed paraffin-embedded BLCA tissue specimen was kindly provided by the Department of Pathology of Nanfang Hospital. The diagnosis of the tumor-node-metastasis (TNM) stages of the BLCA samples depended on the eighth TNM staging system defined by the American Joint Commission on Cancer.

Identification of B cell marker genes

Compared with other public scRNA-seq datasets extracted from GEO, such as GSE135337 [17] and GSE190888 [18], the GSE145137 [19] dataset contained a higher infiltrating ambulance of B cells and was chosen for the B cell markers’ identification. Here, to avoid latent interference, we only used the tumor sample isolated from a patient with primary BLCA (GSM4307111) to conduct the analyses, and a total of 2075 cell samples were included. The Seurat package of R software was utilized to process the scRNA-seq data. Principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE) were conducted to divide the cell samples into various cell clusters. The cell marker genes with identified by the “FindAllMarkers” function of the Seurat package, and the filtering thresholds were set as follows: | logarithmic fold change [logFC]| > 0.5, adjusted P value < 0.05, and pct > 0.25, where pct represented the proportion of a gene with non-zero values in a cell cluster. At last, the SingleR package in R was utilized to perform the cell type annotation.

B cells’ infiltrating profiles identified by computational biology and TIDE algorithm

Different publicly available algorithms, including xCELL, TIMER, quanTIseq, MCPCOUNTER, and EPIC, were employed to estimate the B cells’ infiltration proportion. TIDE algorithm, aiming to predict the immunotherapeutic response based on the transcriptome profiles of cancer samples, was adopted to partly verify the predictive ability of the established model to immunotherapeutic sensitivity [20].

Genomic difference analysis, gene functional annotation, and protein-protein interaction (PPI) network construction

The edgeR package was implemented to identify the differentially expressed genes between paracarinoma and BLCA tissues with |logFC| > 0.5 and false discovery rate (FDR) < 0.05 filtering. The functional annotation of the screened genes was performed through the Metascape (https://metascape.org) with default settings. The PPI network was constructed in the STRING database (https://cn.string-db.org) with a confidence level = 0.4.

Meta-analyses

The meta-analyses were conducted to pool the hazard ratios (HRs) through the meta package in R. As a general rule, if the heterogeneity test was not significant (I² >50% and P value < 0.05), the random effect model would be used to generate the summary effect size of studies; otherwise (I² <50% and P value > 0.05), the fixed effect model would be employed to combine the effect values.

Risk model construction

The genes con-determined by the B cell marker genes’ identification and the genomic divergence analysis between adjacent normal and tumor tissues were chosen for further study. As previously described [21], a gene-pair strategy was employed to construct numerous gene pairs to render the predictive model applicable in different gene expression detection platforms. Briefly, “gene A|gene B” represented a gene pair, and this combination would be regarded as 1 if the expression of gene A is higher than that of gene B; conversely, it would be considered as 0. LASSO regression with 10-fold cross-validation and random forest with nsplit at 10 at the variable hunting was used to identify the significant gene pairs associated with the overall survival (OS) through the glmnet and the randomforestSRC packages, respectively. In the random forest analysis, the VIMP of each variable was calculated, which was positively associated with the variable importance, and the pair with VIMP > 0.001 was chosen. At last, the predictive model was constructed via the multivariate Cox regression with stepwise, and the risk score, named B cell-related score (BCRS), of each sample was calculated. The BCRS was constructed according to the following formula: $\text{BCRS}={{\displaystyle \sum }}_{\text{i}=1}^{\text{n}}{\text{Coeff}}_{\text{i}}\ast {(\text{gene} \text{pair})}_{\text{i}}$. The optimal cut-off value to divide the patients in the TCGA-BLCA cohort into a high- or low-BCRS subgroup was detected by the X-tile, and the cut-off was then applied in all the cohorts enrolled in this study. The survival package in R was used to conduct the Kaplan-Meier survival analyses with log-rank tests. The nomogram to visualize the predictive model was drawn by the rms package.

Pan-cancer analyses

The transcriptome profiles and their corresponding survival information of the other 32 cancer types in the TCGA database were obtained from the UCSC Xena website (https://xenabrowser.net/datapages/). The survival information included disease-free interval (DFI), disease-specific survival (DFS), OS, and progression-free interval (PFI), and the cases with less than 30 days’ follow-up were excluded.

Gene set enrichment analysis (GSEA)

GSEA was performed based on the GSEA software (version 4.3.2), which was directly downloaded from the GSEA official website (http://www.gsea-msigdb.org/gsea/). The reference gene sets about the biological processes of B cells were obtained from the Molecular Signature Database (http://www.gsea-msigdb.org/gsea/msigdb/). The terms with Nominal P < 0.05 and FDR < 0.25 were considered to be statistically significant.

Cell culture

The human normal urethral epithelial cell line, SVHUC1, and two human BLCA cell lines, including T24, and UMUC3, were purchased from the Type Culture Collection of the Chinese Academy of Sciences (Shanghai, China). SVHUC1 cells were maintained in Ham's F-12K Media (Invitrogen, USA). T24 and UMUC3 cells were incubated in Dulbecco’s Modified Eagle Medium (Invitrogen, USA). Those media were supplemented with 10% FBS and 1% antibiotics, and all of the cells were cultured in a humidified atmosphere with 5% CO2 at 37° C.

Real-time quantitative PCR (RT-qPCR)

Trizol reagent (ThermoFisher Scientific, Germany) was used to extract the total RNA of the BLCA and paracarinoma tissue. PrimeScript RT Reagent Kit (Takara, China) and SYBR Premix ExTaq kit (Takara, China) were utilized for cDNA synthesis and amplification. GAPDH was chosen as the internal reference gene, and the detected expression values were normalized via the 2-ΔΔCt based on the ABI Prism 7000 system (Applied Biosystems, USA). The primer sequence used is shown in Supplementary Table 2.

Immunohistochemical staining (IHC) and B cell counting

IHC was conducted in the BLCA specimen using the monoclonal antibody against CD20 (dilution: 1:100, ABclonal, China). The B cells in the slides were distinguished from other cells by staining of CD20, as reported by previous studies [22, 23]. The detailed processes of IHC in the BLCA tissue slides were well described in our earlier study [24]. Each sample was assessed under a microscope at x200 magnification, and the average B cell number was evaluated based on 5 areas with the highest number of CD20+ cells.

Statistical analyses

The statistical analyses of the whole study were conducted via the R software (version 4.0.3) and GraphPad Prism (version 8.0.1). Unless otherwise specified, Wilcoxon signed-rank test was adopted for the continuous variables in two groups, and the Kruskal-Wallis test would be utilized if the number of groups ≥ 3. Pearson’s chi-square test and Fisher’s exact test were performed to measure the divergence for the categorical variables. The Spearman-correlation calculation was performed using the “cor.test” function in R. Welch’s corrected t-test was used to compare the results from the RT-qPCR and B cell counting experiments. The time-dependent receiver operating curve (ROC) analyses were performed and areas under curves (AUCs) were calculated by the timeROC package. In this study, P value < 0.05 was considered significant.

Data availability

The data used to develop the predictive model was downloaded from the TCGA database (https://portal.gdc.cancer.gov/). The validation datasets, including GSE13507, GSE31684, GSE32894, GSE111636, and GSE145137, were obtained in the GEO website (https://ncbi.nlm.nih.gov/geo/). The IMvigor210 dataset was extracted in the IMvigor210CoreBiologies package in R, which could be downloaded from the official website (https://www.r-project.org/). The R code would be provided by the corresponding author upon reasonable request.

Computational biology indicated B cells’ infiltrating was associated with prognosis and immunotherapeutic outcomes in BLCA

The workflow of the whole study is graphically presented in Figure 1. As described above, B cells play a dual role in tumorigenesis. On the one hand, B cells can kill the tumor cells by secreting antibodies and inducing a series of anti-tumor immune processes; on the other hand, B cells also have a pro-tumor effect, which was mainly contributed by the Breg cells and circulating immune complexes secreted by plasma cells [6] (Figure 2A). Besides, despite the fact that the previous studies have disclosed the prognosis value of B cells in 19 cancers, it remains unclear in BLCA [6]. Hence, it is demanded and meaningful to clarify the predictive ability of B cells to clinical outcomes in BLCA at the initial stage of this research. The B cells’ infiltrating profiles were evaluated by multiple publicly available algorithms, including xCELL, TIMER, quanTIseq, MCPCOUNTER, and EPIC, where CIBERSORT-ABS and CIBERSORT were not adopted because they only analyzed the B cells’ subtypes. The infiltration levels in each BLCA sample in the TCGA-BLCA cohort are shown in Supplementary Table 3. The meta-analysis indicated that the subjects with high B cells’ infiltration exhibited a favorable OS rate (pooled HR = 0.84, 95% confidence interval [CI] = 0.72-0.98, Figure 2B). After evaluating the immunotherapeutic response through TIDE, we observed that B cell infiltration was a promising predictor for immunotherapeutic sensitivity (Figure 2C). It should be emphasized that the conclusion was drawn by in silico analysis, though it has been verified in many other cancers [6].

Figure 1. The workflow of the present study.

Figure 2. B cells’ prognosis value in BLCA. (A) The dual roles of B cells in tumor immunity. (B) The B cells’ infiltration level was a predictor for OS in the TCGA-BLCA cohort. (C) The predictive ability of B cells’ infiltrating proportion to immunotherapeutic sensitivity. (D) B cell number in each BLCA specimen was quantified using CD20 staining. (E) The B cell counts of the patients form the NH cohort. (F) The survival statuses and follow-up duration of the patients in the NH cohort. (G) High infiltration level of B cells heralded favorable prognosis in the NH cohort. (H–J) The comparison of the predictive ability of B cell number and routine clinicopathological traits to 6- (H), 12- (I), and 12-months’ OS (J). BLCA, bladder cancer; OS, overall survival.

High infiltrating of B cells heralded favorable prognosis in the NH cohort

To better confirm the prognostic value of B cell infiltrating levels, we collected the tumor samples of 19 BLCA patients from the Nanfang Hospital of Southern Medical University, and then performed the CD20 staining to quantify the B cell number in each specimen (Figure 2D). The clinicopathological features of these patients were extracted from the local hospital’s electronic medical record system (Table 1), and the average B cell number of each BLCA sample is displayed in Figure 2E. Additionally, we also recorded the follow-up data and survival statuses of these 19 BLCA patients (Figure 2F). According to the optimal cut-off value detected by the X-tile software, 19 BLCA patients were divided into high- and low-B cell count subgroups, and the patients with lower B cell infiltration levels suffered poorer OS (P < 0.05, Figure 2G). Compared with routine clinicopathological parameters, including age, TNM stages, pathological T stages, and pathological N stages, B cell number in the BLCA samples was of a stronger predictive ability to prognosis, no matter for 6- (Figure 2H), 12- (Figure 2I), or 18-months’ OS (Figure 2J).

Construction of BCRS

2075 cell samples extracted from a BLCA patient were mainly divided into 9 cell types, which are displayed in Figure 3A. A sum of 472 cell marker genes of B cells was identified, which were differentially expressed in B cells compared with other cell types (Supplementary Table 4). At the same time, 13012 differentially expressed genes between the paracarinoma and BLCA samples were screened from the TCGA-BLCA cohort (Figure 3B and Supplementary Table 5), 163 of which also acted as B cell markers and were selected for further investigation (Figure 3C). The functional annotation indicated that the 163 genes were primarily involved in the immune and inflammatory responses (Figure 3D), agreeing with the previous findings. The PPI network of these 163 genes is shown in Figure 3E, implying the tight regulatory relationship of these genes.

Figure 3. Identification of B cell-related genes. (A) 2075 cells from the primary BLCA were divided into 9 different cell types. (B) 13012 genes showing transcriptome expression difference between paracarinoma and BLCA samples were identified in the TCGA-BLCA cohort. (C) 163 genes were con-determined by the B cell marker genes’ analyses and genomic difference analyses. (D) The functional annotation of the 163 genes. (E) The protein-protein interaction network of the 163 genes.

Subsequently, the 163 genes were cyclically singly paired, and 4070 gene pairs were constructed. The workflow of the predictive model construction is graphically displayed in Figure 4A. First, LASSO regression identified 16 gene pairs as the significant predictor for OS in the TCGA-BLCA cohort (Figure 4B), and their coefficients were shown in Figure 4C. Meanwhile, in the random forest analysis, 38 gene pairs were determined (Figure 4D and Supplementary Table 6). 6 gene pairs intersected with the LASSO regression results and random forest outcomes were then included in the multivariate Cox regression with stepwise (Figure 4E). Accordingly, the predictive model was constructed. BCRS was calculated as follows: BCRS = 0.548*(ANXA1|BTG1) + 0.294*(ANXA1|MDK) + 0.645*(CALD1|IL32) + 0.461*(EMP1|ARPC5) - 0.334*(ABRACL|PSMD2). A sum of 9 genes was contained in this predictive model, including ANXA1, CALD1, EMP1, ABRACL, BTG1, MDK, IL32, ARPC5, and PSMD2, almost all of which were significantly associated with B cell infiltration proportion (all P < 0.05, Figure 4F), indicating that BCRS might be associated with B cell infiltration in BLCA samples (see below for more analyses). Figure 4G shows the expression levels of these 9 genes in each cell cluster identified by the scRNA-seq analysis of the BLCA sample.

Figure 4. Establishment of the B cell-related gene signature. (A) The workflow showed the process of the predictive model construction. (B) The LASSO regression identified 16 gene-pairs significantly associated with OS in the TCGA-BLCA cohort. (C) The coefficients of the variables in the LASSO regression model. (D) The random forest algorithm identified 38 gene-pairs with VIMP > 0.001. (E) 6 gene-pairs were con-determined by the LASSO and random forest analyses, 5 of which were identified with the multivariate Cox regression with stepwise. (F) The association of the genes in the predictive model with B cells’ infiltration proportion. (G) The expression levels of the genes in the established model in each cell cluster.

BCRS was a robust prognosis predictor for BLCA

The HRs and the 95% confidence intervals (CIs) of the variables in the established Cox proportional hazard model are shown in Figure 5A. To help the researchers and clinicians better understand this model, we established a nomographic chart (Figure 5B), which transformed the complex multiple gene-pair model into a visual graph. Next, we used various public external cohorts to confirm the prognosis value of BCRS. The baseline clinical information of the training and external validation cohorts is presented in Table 1. According to the same cut-off, which equaled 2.75, the subjects from the TCGA-BLCA cohort, GSE13507 cohort, GSE31684 cohort, and GSE32894 cohort were divided into the high- or low-BCRS subgroups. The distribution of the BCRS and the survival statuses of the cases from the TCGA-BLCA cohort, GSE13507 cohort, GSE31684 cohort, and GSE32894 cohort are shown in Figure 5C–5F, respectively. Kaplan-Meier survival analyses displayed the predictive ability of BCRS to OS in the TCGA-BLCA cohort (P < 0.001, Figure 5G), GSE13507 cohort (P < 0.01, Figure 5H), GSE31684 cohort (P > 0.05 Figure 5I), and GSE32894 cohort (P < 0.001, Figure 5J). Besides, more death could be observed in the high-BCRS subjects from these cohorts, as shown in Figure 5K–5N.

Figure 5. The validation of the prognosis value of BCRS in BLCA. (A) The hazard ratios and the 95% confidence intervals of the variables in the predictive model. (B) A nomogram was drawn to visualize the predictive model. (C–F) The distribution of the BCRSs and survival statuses of the BLCA cases from the TCGA-BLCA cohort (C), the GSE13507 cohort (D), the GSE31684 cohort (E), and the GSE32894 cohort (F). The cut-off values to divide all the subjects in these cohorts into low- and high-BCRS subgroup was equal to 2.75, which was detected by X-tile in the training cohort. (G–J) The Kaplan-Meier survival analyses of BCRS in the TCGA-BLCA cohort (G), the GSE13507 cohort (H), the GSE31684 cohort (I), and the GSE32894 cohort (J). (K–N) The association between the survival statuses and BCRS stratification in the TCGA-BLCA cohort (K), the GSE13507 cohort (L), the GSE31684 cohort (M), and the GSE32894 cohort (N). BCRS, B cell-related score.

Besides, 19 BLCA samples collected from the local hospital were used for re-verification. The real-time quantitative PCR experiments were conducted to measure the expression values of the genes included in the model (Figure 6A), and the raw CT levels are shown in Supplementary Table 7. The concrete expression values of the genes and the corresponding clinicopathological information of the patients are supplemented in Supplementary Table 8. Through the RT-qPCR experiments, it was observed that ANXA1, CALD1, EMP1, and PSMD2 were up-regulated in the BLCA samples in comparison to the adjacent normal tissue, while IL32 was down-regulated (all P < 0.05, Supplementary Figure 1A). Additionally, we found that the high expression of ANXA1 was significantly associated with unfavorable prognosis in the NH cohort (P < 0.001, Supplementary Figure 1B). To further clarify the association of ANXA1 with the tumorigenesis of BLCA, cellular experiments were then conducted. Compared with SVHUC1 cells, the T24 (P < 0.05) and UMUC3 (P < 0.01) cells exhibited higher levels of ANXA1 (Supplementary Figure 2).

Figure 6. BCRS was a robust and reliable prognosis predictor in BLCA. (A) The transcriptome expression levels of the model’s genes in the local BLCA samples, which was detected by real-time quantitative PCR. (B) A 0-or-1 matrix was established according to the mRNA expression values of the patients. (C) The BCRSs of the BLCA samples collected from the local hospital. (D) The distribution of survival statuses, follow-up duration, and BCRSs in the NH cohort. According to the same cut-off value (2.75), 19 cases were divided into low- and high-BCRS subgroup. (E) The patients with high BCRS exhibited unfavorable OS rates. (F) The comparison of BCRS and routine clinicopathological features to 6-. 12-, and 18-months’ OS. (G, H) Meta-analyses indicated that BCRS was a significant prognosis predictor, both using the continuous variables (G) and binary variables (H) for analyses.

Using the gene-pair strategy and mRNA expression levels detected by the RT-qPCR experiments, a 0-or-1 matrix of the 5 gene pairs in the predictive model was constructed (Figure 6B). According to the formula mentioned above, the BCRS of each patient was calculated (Figure 6C). The distribution of the survival status, follow-up duration, and BCRSs of the patients is displayed in Figure 6D. The patients in the NH cohort were classified into low- and high-BCRS subgroups based on the cut-off value determined above (2.75), and the cases with high BCRS levels exhibited poorer OS (P < 0.05, Figure 6E). We also compared the prognosis predictive performance of BCRS and routine clinicopathological features to 6-, 12-, and 18-months’ OS, and BCRS showed the highest predictive capability (Figure 6F).

Given that the BCRS was not a significant predictor in the GSE31684 cohort (P > 0.05, Figure 5I), we performed meta-analyses to combine the effects, and the result indicated that BCRS was a significant biomarker for the OS in BLCA, no matter using the continuous variables (pooled HR = 1.37, 95%CI = 1.29-1.46, Figure 6G) and the binary variables (pooled HR = 2.79, 95%CI = 2.22-3.49, Figure 6H).

BCRS was superior to routine clinicopathological features in OS prediction

The predictive ability of BCRS was compared with the routine clinicopathological parameters, including age, gender, tumor grade, TNM stage, pathological T stage, pathological N stage, and M stage. To ensure comparability, we transformed all the continuous parameters into binary variables. The optimal cut-off for age was determined by X-tile. The univariate and multivariate Cox regression indicated that BCRS was an independent prognosis predictor both in the TCGA-BLCA cohort and the NH cohort (all P < 0.05, Supplementary Table 9). Besides, compared with the common clinical features, BCRS exhibited a stronger predictive ability to 1- (Supplementary Figure 3A), 2- (Supplementary Figure 3B), 3- (Supplementary Figure 3C), 4- (Supplementary Figure 3D), and 5-year’s (Supplementary Figure 3E) OS in the TCGA-BLCA cohort.

BCRS was a promising tool to evaluate the immunotherapeutic response in BLCA

The association between BCRS and B cells’ biological processes was detected based on the 326 B cell samples extracted from the scRNAseqDB (https://bioinfo.uth.edu/scrnaseqdb/). Subsequently, GSEA was performed after the calculation of the BCRS of each B cell sample. The detailed processes are displayed in Figure 7A. Since all the BCRSs of these B cells were less than 2.75, we chose the median BCRS level (equaled 1) to classify the B cell samples into low- and high-BCRS subgroups. It was observed that the negative regulation of B cells’ activation (Nominal P < 0.05, FDR < 0.25, Figure 7B) and proliferation (Nominal P > 0.05, FDR < 0.25, Figure 7C) was enriched in the high-BCRS B cells, while the positive regulation of B cells’ activation (Nominal P < 0.05, FDR < 0.25, Figure 7D) and proliferation (Nominal P < 0.05, FDR < 0.25, Figure 7E) was significant in the low-BCRS B cells, suggesting that BCRS was negatively associated with B cells’ activation and proliferation. In the BLCA patients, the cases with high BCRS exhibited low B cell infiltration levels, both in the TCGA-BLCA cohort (Figure 7F) and the NH cohort (P < 0.001, Figure 7G).

Figure 7. BCRS was associated with B cells’ activation and proliferation, and served as a potential predictor to immunotherapeutic response. (A) The processes of the GSEA in the 326 B cell samples. (B, C) The negative regulation of B cells’ activation (B) and proliferation (C) was enriched in the high-BCRS B cells. (D, E) The positive regulation of B cells’ activation (D) and proliferation (E) was enriched in the low-BCRS B cells. (F, G) BCRS was negatively associated with B cells’ infiltration levels, both in the TCGA-BLCA cohort (F) and the NH cohort (G). (H) The levels of BCRS between immunotherapeutic responders and non-responders from the TCGA-BLCA cohort, IMvigor210 cohort, and GSE111636 cohort. (I) The predictive ability of BCRS to immunotherapeutic response in the TCGA-BLCA cohort, IMvigor210 cohort, and GSE111636 cohort. (J) High levels of BCRS heralded poor prognosis in the IMvigor210 cohort.

Next, the immunotherapeutic responses of the cases from the TCGA-BLCA cohort were evaluated by the TIDE algorithm. The subjects with low BCRSs were more likely to respond to the immunotherapy through the Wilcoxon signed-rank test (P < 0.01, Figure 7H). Two BLCA cohorts receiving anti-PD1/PDL1 treatment (pembrolizumab for GSE111636 and atezolizumab for IMvigor210) were also utilized for validation in real-world situations. The cases with high BCRSs were less likely to benefit from the atezolizumab treatment (P < 0.01, Figure 7H) in the IMvigor210 cohort. The levels of BCRS in the patients showing response to pembrolizumab were lower than those in the cases exhibiting non-response, but the difference was not statistically significant (P > 0.05 Figure 7H), partly due to the limited sample size. The ROC analyses indicated the predictive ability of BCRS to immunotherapeutic effectiveness in the GSE111636 cohort, IMvigor210cohort, and the TCGA-BLCA cohort (Figure 7I). Besides, in the 195 BLCA patients from the IMvigor210 cohort, the subjects with high BCRS suffered poor prognoses (P < 0.05, Figure 7J). Totally, BCRS was a promising clinical to distinguish the BLCA patients who were more likely to benefit from the immunotherapy.

Pan-cancer analyses of BCRS

The prognosis value of BCRS at a pan-caner level was also detected. Besides the BLCA, 32 other cancers from the TCGA database were used. Their abbreviations, full terms, and sample size are displayed in Supplementary Table 10. Based on the established formula, the BCRS of each sample was calculated (Supplementary Table 11). It could be observed that BCRSs in different cancers were of high heterogeneity (P < 0.001, Supplementary Figure 4A). The predictive performance of BCRS to DFI, DSS, OS, and PFI in each cancer was then evaluated through univariate Cox regression, and we found that BCRS was a significant prognosis predictor in 21 of 33 cancers (Supplementary Figure 4B and Table 2). At a pan-cancer level, BCRS was negatively associated with B cells’ infiltration proportion in 21 cancer types, and only 6 cancer types showed a positive association (Supplementary Figure 5A). Since tumor mutation burden and microsatellite instability have been widely accepted as biomarkers of immunotherapy, we supplemented the correlation of BCRS to tumor mutation burden (Supplementary Figure 5B) and microsatellite instability (Supplementary Figure 5C) as a reference. Totally, these results implied the tremendous potential of BCRS’s application in other cancer types.