Advertisement
Urothelial Cancer| Volume 51, P78-88, May 2023

Download started.

Ok

Single-nucleus and Spatially Resolved Intratumor Subtype Heterogeneity in Bladder Cancer

Open AccessPublished:April 07, 2023DOI:https://doi.org/10.1016/j.euros.2023.03.006

      Abstract

      Background

      Current bulk transcriptomic classification systems for bladder cancer do not consider the level of intratumor subtype heterogeneity.

      Objective

      To investigate the extent and possible clinical impact of intratumor subtype heterogeneity across early and more advanced stages of bladder cancer.

      Design, setting, and participants

      We performed single-nucleus RNA sequencing (RNA-seq) of 48 bladder tumors and additional spatial transcriptomics for four of these tumors. Total bulk RNA-seq and spatial proteomics data were available from the same tumors for comparison, along with detailed clinical follow-up of the patients.

      Outcome measurements and statistical analysis

      The primary outcome was progression-free survival for non–muscle-invasive bladder cancer. Cox regression analysis, log-rank tests, Wilcoxon rank-sum tests, Spearman correlation, and Pearson correlation were used for statistical analysis.

      Results and limitations

      We found that the tumors exhibited varying levels of intratumor subtype heterogeneity and that the level of subtype heterogeneity can be estimated from both single-nucleus and bulk RNA-seq data, with high concordance between the two. We found that a higher class 2a weight estimated from bulk RNA-seq data is associated with worse outcome for patients with molecular high-risk class 2a tumors. The sparsity of the data generated using the DroNc-seq sequencing protocol is a limitation.

      Conclusions

      Our results indicate that discrete subtype assignments from bulk RNA-seq data may lack biological granularity and that continuous class scores may improve clinical risk stratification of patients with bladder cancer.

      Patient summary

      We found that several molecular subtypes can exist within a single bladder tumor and that continuous subtype scores can be used to identify a subgroup of patients with poor outcomes. Use of these subtype scores may improve risk stratification for patients with bladder cancer, which can help in making decisions on treatment.

      Keywords

      1. Introduction

      Transcriptomic subtyping of bladder cancer (BC) using bulk gene expression profiling has been in focus for two decades [
      • Dyrskjøt L.
      • Thykjaer T.
      • Kruhøffer M.
      • et al.
      Identifying distinct classes of bladder carcinoma using microarrays.
      ]. However, while current bulk classification systems assign a single subtype to each tumor [
      • Kamoun A.
      • de Reyniès A.
      • Allory Y.
      • et al.
      A consensus molecular classification of muscle-invasive bladder cancer.
      ,
      • Lindskrog S.V.
      • Prip F.
      • Lamy P.
      • et al.
      An integrated multi-omics analysis identifies prognostic molecular subtypes of non-muscle-invasive bladder cancer.
      ], expression profiling of multiple tumor regions and immunohistochemistry analyses have shown that different subtypes can coexist within a tumor [
      • Thomsen M.B.H.
      • Nordentoft I.
      • Lamy P.
      • et al.
      Comprehensive multiregional analysis of molecular heterogeneity in bladder cancer.
      ,
      • Sjödahl G.
      • Eriksson P.
      • Lövgren K.
      • et al.
      Discordant molecular subtype classification in the basal-squamous subtype of bladder tumors and matched lymph-node metastases.
      ].
      Single-cell technologies provide an opportunity to study tumor ecosystems and heterogeneity at single-cell resolution. Single-cell RNA sequencing (scRNA-seq) requires cells isolated from fresh tissue, whereas single-nucleus RNA sequencing (snRNA-seq) makes it possible to investigate frozen tumor tissue, facilitating biobank studies of clinically well-annotated patients. Studies of heterogeneity within the epithelial compartment of BC are emerging [
      • Sfakianos J.P.
      • Daza J.
      • Hu Y.
      • et al.
      Epithelial plasticity can generate multi-lineage phenotypes in human and murine bladder cancers.
      ,
      • Lai H.
      • Cheng X.
      • Liu Q.
      • et al.
      Single-cell RNA sequencing reveals the epithelial cell heterogeneity and invasive subpopulation in human bladder cancer.
      ,
      • Wang H.
      • Mei Y.
      • Luo C.
      • et al.
      Single-cell analyses reveal mechanisms of cancer stem cell maintenance and epithelial-mesenchymal transition in recurrent bladder cancer.
      ,
      • Gouin 3rd, K.H.
      • Ing N.
      • Plummer J.T.
      • et al.
      An N-Cadherin 2 expressing epithelial cell subpopulation predicts response to surgery, chemotherapy and immunotherapy in bladder cancer.
      ]; however, an investigation of the extent and possible clinical impact of intratumor heterogeneity at single-cell resolution in larger cohorts of both early and advanced BC stages is still lacking.
      We performed snRNA-seq analysis of 59 052 nuclei from 48 bladder tumors to study tumor heterogeneity at single-cell resolution. We demonstrate that the tumors exhibit varying levels of intratumor subtype heterogeneity, which may affect clinical outcomes.

      2. Materials and methods

      The Supplementary material provides full details of the analyses performed.

      2.1 Tumor samples

      A total of 48 bladder tumor samples (10 Ta, 13 T1, and 25 T2–4) were obtained via transurethral resection of the bladder or radical cystectomy (Supplementary Table 1). One sample was processed directly from surgery (patient 4059) and 47 samples were embedded in O.C.T., frozen in liquid nitrogen and stored at −80°C.

      2.2 snRNA-seq using an optimized DroNc-seq protocol

      We performed snRNA-seq of bladder tumors using an optimized DroNc-seq protocol previously described [
      • Schmøkel S.S.
      • Nordentoft I.K.
      • Lindskrog S.V.
      • et al.
      Improved protocol for single nucleus RNA-sequencing of frozen human bladder tumor biopsies.
      ]. Seurat v3.2.0 [
      • Stuart T.
      • Butler A.
      • Hoffman P.
      • et al.
      Comprehensive integration of single-cell data.
      ] was used to process and analyze the data. All epithelial nuclei were classified according to the UROMOL2021 classes [
      • Lindskrog S.V.
      • Prip F.
      • Lamy P.
      • et al.
      An integrated multi-omics analysis identifies prognostic molecular subtypes of non-muscle-invasive bladder cancer.
      ] or the consensus classes for muscle-invasive BC (MIBC) [
      • Kamoun A.
      • de Reyniès A.
      • Allory Y.
      • et al.
      A consensus molecular classification of muscle-invasive bladder cancer.
      ].

      2.3 RNA-seq

      Bulk RNA-seq data were available for 44 tumors. We used the R package Weighted In-Silico Pathology (WISP) v2.3 [

      WISP. Weighted In Silico Pathology: a novel approach to assess intra-tumoral heterogeneity. https://cit-bioinfo.github.io/WISP/.

      ] to deconvolute intratumor heterogeneity from bulk transcriptomic profiles.

      2.4 Spatial transcriptomics

      Four samples were analyzed using 10x Visium Spatial (10x Genomics, Pleasanton, CA, USA). Sections of tissue were cut for methanol fixation, hematoxylin and eosin staining, and imaging, followed by Visium Spatial gene expression library construction according to the manufacturer’s instructions (CG000160 and CG000239). Seurat v4.2.0 and STutility v1.1.1 [
      • Bergenstråhle J.
      • Larsson L.
      • Lundeberg J.
      Seamless integration of image and molecular analysis for spatial transcriptomics workflows.
      ] were used for subsequent data filtering, normalization, and visualization.

      2.5 Statistical analysis

      Progression-free survival (PFS) was assessed for patients with non–muscle-invasive BC (NMIBC) and was defined as the time to progression to MIBC. Survival analysis was performed using the Kaplan-Meier method and two-sided log-rank tests were used to compare survival curves (R packages survival v3.2.13 and survminer v0.4.9). Univariate and multivariable Cox regression analyses of PFS were performed using clinical and molecular predictors. Association between the fraction of lymphocytes in tumors collected before bacillus Calmette-Guérin (BCG) treatment and progression to MIBC was assessed using the Wilcoxon rank-sum test. All analyses were performed with R v3.6.1 (R Foundation for Statistical Computing, Vienna, Austria), except for the SCENIC analysis, which was performed with Python v3.10.5 (Python Software Foundation, Delaware, USA).

      3. Results

      3.1 snRNA-seq of human bladder tumors: data processing and method comparison

      In order to conduct retrospective analyses on clinically well-annotated samples, we performed snRNA-seq of tumors from 48 patients with BC (Supplementary Table 1) using an optimized DroNc-seq protocol [
      • Schmøkel S.S.
      • Nordentoft I.K.
      • Lindskrog S.V.
      • et al.
      Improved protocol for single nucleus RNA-sequencing of frozen human bladder tumor biopsies.
      ]. We obtained data for 117 653 nuclei after initial barcode processing (Fig. 1A), of which 59 052 nuclei passed downstream quality control filtering (Fig. 1B).
      Figure thumbnail gr1
      Fig. 1Single-nucleus RNA sequencing of human bladder tumors: data processing and method comparison. (A) Distribution of the number of genes (left) and counts (right) detected per nucleus for data sets with different Hamming distances used to correct barcodes. (B) Number of genes (left) and counts (right) per nucleus for the full DroNc-seq data set after quality control filtering (n = 59 052). (C) Number of genes and counts per nucleus for three samples analyzed using the 10x Chromium platform (n = 5758) and the DroNc-seq platform (n = 4431). (D) UMAP visualization of 5010 nuclei from three samples analyzed using 10x Chromium, colored by patient origin (left) and cell type (right). (E) UMAP visualization of 3903 nuclei from three samples analyzed using DroNc-seq, colored by patient origin (left) and cell type (right). (F) Relative proportions of cell types within three samples analyzed on two platforms: 10x Chromium (2004, n = 2569; 2514, n = 895; 4735, n = 1546) and DroNc-seq (2004, n = 1525; 2514, n = 1794; 4735, n = 584). UMAP = uniform manifold approximation and projection.
      To investigate the quality of data generated using DroNc-seq, three samples were subsequently analyzed using 10x Chromium. After quality control filtering, 4431 nuclei (average of 568 expressed genes per nucleus) remained when using DroNc-seq, whereas 5758 nuclei (average of 1649 expressed genes per nucleus) remained when using 10x Chromium (Fig. 1C). For genes with nonzero expression in both the DroNc-seq and 10x Chromium data sets, we observed 193 and 1054 counts per gene on average, respectively. Thus, 10x Chromium generated more data per nucleus in comparison to DroNc-seq, but at a higher cost. We assigned all nuclei to one of five major cell populations (epithelial, fibroblast, endothelial, lymphoid, or myeloid) using data from a previous study of snRNA-seq of human bladder tumors as the reference [
      • Gouin 3rd, K.H.
      • Ing N.
      • Plummer J.T.
      • et al.
      An N-Cadherin 2 expressing epithelial cell subpopulation predicts response to surgery, chemotherapy and immunotherapy in bladder cancer.
      ]. Uniform manifold approximation and projection (UMAP) visualization of the 10x Chromium data revealed three major epithelial clusters separated by patient origin, as well as smaller clusters consisting of nonepithelial cells (Fig. 1D). By contrast, nonepithelial cells clustered together with the patient-specific epithelial cells for DroNc-seq data, probably because of the lower number of genes detected per nucleus (Fig. 1E). To confirm cell-type annotations for the more sparse DroNc-seq data, we compared the relative proportions of cell-type annotations within each tumor obtained via the two platforms. We found a similar proportion of nonepithelial cells, confirming that the epithelial compartment constitutes the bulk of these tumors when using either platform (Fig. 1F). We did observe differences in the fractions of nonepithelial cell types between the platforms, which may be caused by tumor heterogeneity (not the same cells analyzed) and by differences in the number of genes detected per nucleus (Supplementary Fig. 1A). Although comparison of the methods revealed better performance with 10x Chromium, DroNc-seq allowed us to analyze a larger cohort of bladder tumors at single-cell resolution because of the lower cost of the analysis.

      3.2 Exploration of the tumor microenvironment

      The DroNc-seq data set consisted of 48 bladder tumors spanning the BC disease spectrum (10 Ta, 13 T1, and 25 T2–4 tumors). In total, 48 567 nuclei were annotated to a major cell type, and the various cell populations showed higher expression of corresponding marker genes (Fig. 2A). Overall, the 48 tumors were predicted to consist of 91% epithelial (cancer) cells, 5% fibroblasts, 3% immune cells, and 1% endothelial cells. Graph-based clustering of all tumors was mainly driven by patient origin, indicating a high level of intertumor heterogeneity (Fig. 2B). The fraction of nonepithelial cells varied between tumors and increased with tumor stage, as expected (Fig. 2C). Other data layers, including multiplex immunofluorescence (mIF) and bulk RNA-seq, have previously been generated from the same tumor samples, facilitating comparison of the estimated immune-cell percentage across these layers [
      • Taber A.
      • Christensen E.
      • Lamy P.
      • et al.
      Molecular correlates of cisplatin-based chemotherapy response in muscle invasive bladder cancer by integrated multi-omics analysis.
      ,
      • Taber A.
      • Prip F.
      • Lamy P.
      • et al.
      Immune contexture and differentiation features predict outcome in bladder cancer.
      ,
      • Strandgaard T.
      • Lindskrog S.V.
      • Nordentoft I.
      • et al.
      Elevated T-cell exhaustion and urinary tumor DNA levels are associated with bacillus Calmette-Guérin failure in patients with non-muscle-invasive bladder cancer.
      ]. We observed that the majority of tumors with a higher number of immune cells estimated from DroNc-seq data also had a higher percentage of immune cells estimated from mIF data and a relatively higher bulk RNA-seq immune score (Fig. 2D,E).
      Figure thumbnail gr2
      Fig. 2Exploration of the tumor microenvironment of human bladder tumors. (A) Violin plots showing the mean expression of epithelial, fibroblast, endothelial, lymphoid, and myeloid gene markers in the different cell types. Cell-type markers were identified from the reference data set using SingleR. (B) UMAP visualization of 48 567 nuclei from 48 tumors analyzed using the DroNc-seq platform, colored by patient origin, tumor stage, tumor grade, and cell type. (C) Relative proportions of cell types in each sample, separated by tumor stage. Samples are sorted by increasing number of nuclei in each tumor stage. (D) Correlation between the percentage of immune cells from DroNc-seq data and the immune cell percentage from multiplex immunofluorescence (mIF) data for 13 samples. Samples are colored by tumor stage. Spearman correlation was used to determine the correlation coefficient R and p value. (E) Correlation between the percentage of immune cells from DroNc-seq data and the estimated immune score from bulk RNA-seq data for 44 samples. Samples are colored by tumor stage. Spearman correlation was used to determine the correlation coefficient R and p value. UMAP = uniform manifold approximation and projection.
      The strength of the snRNA-seq method is the ability to analyze frozen tumors from patients with known clinical outcomes. Among the patients included in this study, 15 received at least five BCG instillations. Although the data are based on just a few samples and a low number of lymphocytes identified, we found that a lower fraction of lymphocytes in pre-BCG tumors was significantly associated with progression to MIBC (p = 0.018). No significant associations between specific cell populations and BCG or chemotherapy response were observed.

      3.3 Intratumor subtype heterogeneity assessed at single-cell resolution

      We focused our analysis on the epithelial (cancer cell) compartment for each tumor individually. Bulk RNA-seq data were available for 44 of the tumors, and the tumors were classified according to the consensus classes for MIBC [
      • Kamoun A.
      • de Reyniès A.
      • Allory Y.
      • et al.
      A consensus molecular classification of muscle-invasive bladder cancer.
      ] or the UROMOL2021 classes for NMIBC [
      • Lindskrog S.V.
      • Prip F.
      • Lamy P.
      • et al.
      An integrated multi-omics analysis identifies prognostic molecular subtypes of non-muscle-invasive bladder cancer.
      ], depending on tumor stage. To assign a cancer cell phenotype to each epithelial cell and explore intratumor subtype heterogeneity at the single-nucleus level, we used the bulk classification systems to classify all epithelial nuclei from tumors having >100 nuclei and at least 40% of the respective classifier genes expressed (NMIBC: 29 484 nuclei for 17 tumors; MIBC: 9523 nuclei from ten tumors). Tumors exhibited varying levels of intratumor subtype heterogeneity, with the dominating class at the single-nucleus level constituting 36–97% of all nuclei (mean 69%; Fig. 3A,B), indicating that a discrete subtype assignment from bulk RNA-seq data may lack biological granularity. The dominating class at single-nucleus level was consistent with the transcriptomic class from bulk RNA-seq in 79% (11/14) of NMIBC tumors (not considering class 2b tumors) and 67% (6/9) of MIBC tumors. Pseudo-bulk classifications (generated from snRNA-seq data) are shown in Supplementary Table 1.
      Figure thumbnail gr3
      Fig. 3Intratumor subtype heterogeneity accessed at single-cell resolution. (A) Relative proportion of single-nucleus subtype classifications of epithelial cells in NMIBC tumors separated by the subtype from bulk RNA sequencing. Classification according to the UROMOL2021 system. Samples are sorted by increasing separation levels from the bulk classification. (B) Relative proportion of single-nucleus subtype classifications of epithelial cells in MIBC tumors separated by the subtype from bulk RNA sequencing. Classification according to the consensus MIBC system. Samples are sorted by increasing separation levels from the bulk classification. The tumor from patient 4735 was not classified, as no bulk RNA sequencing was available. (C) Immunohistochemical staining for CK5/6 and GATA3 of three areas (tissue microarray cores) from the tumor from patient 1468. (D–G) UMAP visualization of selected tumors colored by single-nucleus classifications (top) and dot plots showing the mean expression of subtype-specific gene signatures in each subtype population (bottom). The numbers of nuclei within each class were: (D) patient 4735; LumP, 71; LumU, 9; LumNS, 13; Ba/Sq, 1; and NE-like, 465; (E) patient 2704: LumP, 163; LumU, 234; LumNS, 21; Ba/Sq, 14; stroma-rich, 11; and NE-like, 564; (F) patient 3096: LumP, 139; LumU, 44; LumNS, 67; Ba/Sq, 171; stroma-rich, 24; and NE-like, 32; and (G) patient 2004: class 1, 7; class 2a, 435; class 2b, 24; and class 3, 955. Ba/Sq = basal/squamous; LumNS = luminal nonspecified; LumP = luminal papillary; LumU = luminal unstable; MIBC = muscle-invasive bladder cancer; NE = neuroendocrine; NMIBC = non–muscle-invasive bladder cancer; UMAP = uniform manifold approximation and projection
      We further investigated the intratumor heterogeneity using immunohistochemical (IHC) staining of the MIBC tumor from patient 1468. According to the single-nucleus analysis, this tumor was mainly of luminal origin but it also consisted of a minor fraction of nuclei classified as basal/squamous (Ba/Sq; Fig. 3B). Interestingly, IHC staining of three tissue microarray cores revealed areas with malignant cells positive for the basal marker cytokeratin 5/6 (CK5/6), while the majority of malignant cells in all cores were GATA3-positive, confirming the presence of both luminal and basal features simultaneously (Fig. 3C).
      Next, we performed graph-based clustering of the tumors individually and found cases for which the transcriptomic classes were a major driver of the clustering (Fig. 3D–G; Supplementary Fig. 1B,C). A histological neuroendocrine (NE) tumor (patient 4735) consisted of two major populations: nuclei classified as NE-like and nuclei classified as luminal papillary (LumP). This subtype-dependent division was consistent with the expression of marker genes for NE tumors (Fig. 3D). In agreement with this, a tumor classified as NE-like from the bulk RNA-seq (patient 2704) showed a major cluster of mixed nuclei classified as either NE-like or luminal unstable (LumU) and a minor cluster of LumP classifications (Fig. 3E). Tumors from both patient 3096 and patient 2004 consisted of a major cluster in accordance with the transcriptomic class from bulk RNA-seq (Ba/Sq and class 3, respectively), along with other minor subtype clusters, further emphasizing the presence of intratumor subtype heterogeneity beyond discrete bulk classifications (Fig. 3F,G).
      Next, we applied SCENIC [
      • Van de Sande B.
      • Flerin C.
      • Davie K.
      • et al.
      A scalable SCENIC workflow for single-cell gene regulatory network analysis.
      ], a method for reconstructing regulons and inferring regulon activity scores for each nuclei, using all nuclei with >500 genes and tumors with at least 50 nuclei (21 583 nuclei from 39 tumors). We explored regulon activities associated with the single-nucleus subtype classifications for NMIBC and MIBC tumors to further support the classifications. In agreement with findings from the UROMOL2021 study [
      • Lindskrog S.V.
      • Prip F.
      • Lamy P.
      • et al.
      An integrated multi-omics analysis identifies prognostic molecular subtypes of non-muscle-invasive bladder cancer.
      ], class 3 nuclei showed significantly higher activity of the TP63, KDM5B, and RARG regulons (Supplementary Fig. 2A). For MIBC, luminal tumors showed higher activity of, for example, the GATA3 and FOXA1 regulons, and stroma-rich tumors showed higher activity of the PGR and STAT3 regulons, as previously observed (Supplementary Fig. 2B) [
      • Kamoun A.
      • de Reyniès A.
      • Allory Y.
      • et al.
      A consensus molecular classification of muscle-invasive bladder cancer.
      ]. Furthermore, subtype-specific regulon activities were also associated with single-nucleus subtype classifications when analyzing selected tumors individually (Fig. 4A–D).
      Figure thumbnail gr4
      Fig. 4Regulon activity scores for single nuclei. (A–D) UMAP visualization of selected tumors, colored by single-nucleus subtype classifications (left) and activity scores for selected subtype-specific regulons (middle, right). Dot plots of the mean regulon activity in each subtype population are shown to the right of each regulon-colored UMAP. The numbers of nuclei within each class were: (A) patient 2004: class 1, 2; class 2a, 202; class 2b, 16; and class 3, 515; (B) patient 3096: LumP, 45; LumU, 15; LumNS, 14; Ba/Sq, 72; stroma-rich, 14; and NE-like, 4; (C) patient 2704: LumP, 66; LumU, 88; LumNS, 6; Ba/Sq, 6; and NE-like, 205; and (D) patient 4735: LumP, 35; LumU, 1; LumNS, 2; and NE-like, 220. Ba/Sq = basal/squamous; LumNS = luminal nonspecified; LumP = luminal papillary; LumU = luminal unstable; NE = neuroendocrine; UMAP = uniform manifold approximation and projection.

      3.4 Intratumor subtype heterogeneity assessed via spatial transcriptomics

      To explore the spatial organization of intratumor subtype heterogeneity, we performed spatially resolved whole transcriptomics of two T1 tumors (patients 2919 and 3197) and two T2–4 tumors (patients 4735 and 3143) using Visium Spatial. We classified each spot using the appropriate bulk classification system (Fig. 5A–D). The tumors from patients 2919 and 4735 showed high subtype homogeneity at the single-nucleus level, which was confirmed at the spatial transcriptomics level (Fig. 5A,C). Although identified as class 2a at the bulk level, the tumor from patient 3197 showed a high proportion of low-risk class 3 at the single-nucleus level, and low-risk classifications (classes 1 and 3) were also dominant on spatial transcriptomics (Fig. 5B). Finally, the spots for patient 3143, which was classified as LumU at the bulk level, were mainly classified as luminal nonspecified (LumNS). This could be explained by tumor heterogeneity, as only a smaller tissue section was analyzed using spatial transcriptomics (Fig. 5D). However, a minor fraction of the LumNS classification was observed at the single-nucleus level as well (Fig. 3B). Furthermore, we used the anchor-based integration workflow in Seurat to deconvolute the cellular composition of each spot using the same snRNA-seq reference as previously [
      • Gouin 3rd, K.H.
      • Ing N.
      • Plummer J.T.
      • et al.
      An N-Cadherin 2 expressing epithelial cell subpopulation predicts response to surgery, chemotherapy and immunotherapy in bladder cancer.
      ]. In general, we observed that areas with higher scores for nonepithelial cell compartments were classified as class 2b for NMIBC tumors or stroma-rich for MIBC tumors (Fig. 5A–D).
      Figure thumbnail gr5
      Fig. 5Intratumor subtype heterogeneity assessed via spatial transcriptomics. (A–D) Images of hematoxylin and eosin staining of tissue areas (left) from the four tumors analyzed using 10x Visium Spatial. Spots are colored by subtype classifications (middle) and cellular composition scores (right). Ba/Sq = basal/squamous; LumNS = luminal nonspecified; LumP = luminal papillary; LumU = luminal unstable; NE = neuroendocrine.

      3.5 Delineation of intratumor subtype heterogeneity from bulk transcriptomic profiles

      On the basis of this knowledge of intratumor subtype heterogeneity at the single-nucleus level, we used bulk RNA-seq data to investigate the impact of the presence of high-risk subtypes in tumor subclones that may not be identified from bulk subtype assignments. To approximate intratumor subtype heterogeneity from bulk transcriptomic profiles, we used the deconvolution tool WISP [

      WISP. Weighted In Silico Pathology: a novel approach to assess intra-tumoral heterogeneity. https://cit-bioinfo.github.io/WISP/.

      ] on the UROMOL2021 (n = 505) and The Cancer Genome Atlas (TCGA) cohorts (n = 406). As observed at the single-nucleus level, tumors exhibited varying levels of intratumor subtype heterogeneity (Fig. 6A,B). WISP weights correlated with the proportion of class assignments at the single-nucleus level, especially for class 2a (Fig. 6C; Supplementary Fig. 3A). Finally, we explored the clinical impact of the weight of class 2a, as class 2a is the molecular high-risk subtype of NMIBC tumors. We found that a high class 2a weight (estimated using WISP) was significantly associated with worse PFS when considering tumors with a bulk class 2a (p = 0.0005) and class 2b assignment (p = 0.042; Fig. 6D; Supplementary Fig. 3B). Notably, for patients with bulk class 2a tumors, a high class 2a weight was still associated with worse PFS when adjusted for the European Association of Urology clinical risk categories (hazard ratio 2.97, 95% confidence interval [CI] 1.22–7.23; p = 0.017) and European Organization for Research and Treatment of Cancer risk score (hazard ratio 2.58, 95% CI 1.06–6.29; p = 0.036; Supplementary Table 2). For low-risk class 1 and 3 tumors, a high class 2a weight was not associated with worse PFS, but a significant correlation to recurrence-free survival was observed (p = 0.0005; Supplementary Fig. 3B,C).
      Figure thumbnail gr6
      Fig. 6Delineation of intratumor subtype heterogeneity from bulk transcriptomic profiles. (A) WISP weights for the four UROMOL2021 classes for 505 tumors from the UROMOL2021 cohort separated by bulk subtype assignment. (B) WISP weights for the six consensus classes of muscle-invasive bladder cancer (MIBC) for 406 tumors from The Cancer Genome Atlas (TCGA) cohort separated by bulk subtype assignment. (C) Correlation between bulk WISP weights and DroNc-seq class proportions for each of the UROMOL2021 classes. Samples are colored by bulk subtype assignment. Pearson correlation was used to determine the correlation coefficient R and p value. (D) Kaplan-Meier plot of progression-free survival (PFS) for 136 patients with class 2a tumors, stratified by class 2a weight from WISP analysis (log-rank test). An optimal cutpoint for the class 2a WISP weight according to time to progression was determined using the R package survminer (cutpoint 0.5562). Ba/Sq = basal/squamous; LumNS = luminal nonspecified; LumP = luminal papillary; LumU = luminal unstable; NE = neuroendocrine; HR = hazard ratio; CI = confidence interval; RNA-seq = RNA sequencing.

      4. Discussion

      We performed snRNA-seq analysis of tumors from 48 patients with BC. We demonstrated that tumors from both early and more advanced disease stages exhibit large intratumor subtype heterogeneity and that subtype heterogeneity could be estimated from both single-nucleus and bulk RNA-seq data, with high concordance between the two. Notably, we showed that a high class 2a weight estimated from bulk RNA-seq was associated with worse outcome in patients with molecular high-risk class 2a tumors, indicating that discrete subtype assignments from bulk RNA-seq may lack biological granularity. In line with our observations, a previous scRNA-seq study of three murine and two human bladder tumors provided initial evidence of intratumor subtype heterogeneity [
      • Sfakianos J.P.
      • Daza J.
      • Hu Y.
      • et al.
      Epithelial plasticity can generate multi-lineage phenotypes in human and murine bladder cancers.
      ]. Although validation is needed, this supports the potential use of continuous classification scores instead of discrete class assignments for improved risk stratification of patients. Importantly, intratumor subtype heterogeneity may need to be considered in subtype-guided clinical trials. In fact, the BISCAY trial, a biomarker-driven multiarm adaptive trial, showed no differences in response rate between biomarker-directed combination arms and durvalumab monotherapy.
      Previous studies have used scRNA-seq data to refine the current bulk molecular subtypes for breast and colorectal cancer [
      • Wu S.Z.
      • Al-Eryani G.
      • Roden D.L.
      • et al.
      A single-cell and spatially resolved atlas of human breast cancers.
      ,
      • Joanito I.
      • Wirapati P.
      • Zhao N.
      • et al.
      Single-cell and bulk transcriptome sequencing identifies two epithelial tumor cell states and refines the consensus molecular classification of colorectal cancer.
      ], highlighting how single-cell data can improve the biological granularity of subtyping. Owing to the sparsity of the DroNc-seq data presented here, we did not develop new classification systems compatible with snRNA-seq data and did not refine the bulk classification systems for BC. Instead, we applied the existing bulk classification systems to the snRNA-seq data. Inconsistency between the dominating class at single-nucleus level and the class from bulk RNA-seq may therefore be explained by several levels of heterogeneity (different parts of the tumor go to the various analyses) and methodological differences, including the technical challenges of applying bulk classification systems to single-nucleus data. However, we found good concordance overall between bulk and single-nucleus classifications.
      The ability to analyze frozen tumors with known clinical outcomes gives snRNA-seq a major advantage over scRNA-seq. Although we did not thoroughly explore the minor fractions of nonepithelial tumor compartments, we found that a lower fraction of lymphocytes in pre-BCG tumors from 15 patients treated with minimum of five BCG instillations was significantly associated with progression to MIBC. This aligns with previous findings based on mIF staining, where high cytotoxic T-cell infiltration was associated with a lower risk of progression [
      • Taber A.
      • Prip F.
      • Lamy P.
      • et al.
      Immune contexture and differentiation features predict outcome in bladder cancer.
      ]. We observed that the highly infiltrated subtypes (ie, class 2b for NMIBC and stroma-rich for MIBC) were not captured when looking at single nuclei from cancer cells only, as high infiltration characterizes these subtypes. Therefore, given the clinical impact of infiltration, the level of infiltration should still be considered if bulk classification systems are refined using epithelial cell phenotypes observed from single-cell data.
      Finally, our comparison of methods showed better performance of 10x Chromium, with a notably higher level of data per nuclei in comparison to DroNc-seq. This could be because of the higher capture rate, better error suppression in barcodes, or the use of elastic beads provided by 10x Genomics. However, while 10x Chromium is less time-consuming and a more standardized system, the DroNc-seq system has lower running costs and was therefore selected for this exploratory study to increase the cohort size.

      5. Conclusions

      Our results highlight that intratumor subtype heterogeneity is an important biological feature of human bladder tumors that may have an impact on clinical outcome. Continuous subtype scores may have the potential to refine biological characterization of tumors and clinical stratification of patients.
      Author contributions: Lars Dyrskjøt had full access to all the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis.
      Study concept and design: Lindskrog, Schmøkel, Nordentoft, Dyrskjøt.
      Acquisition of data: Schmøkel, Nordentoft, Dyrskjøt.
      Analysis and interpretation of data: Lindskrog, Schmøkel, Nordentoft, Lamy, Knudsen, Prip, Strandgaard, Dyrskjøt.
      Drafting of the manuscript: Lindskrog, Schmøkel, Dyrskjøt.
      Critical revision of the manuscript for important intellectual content: All authors.
      Statistical analysis: Lindskrog, Schmøkel, Lamy.
      Obtaining funding: Dyrskjøt.
      Administrative, technical, or material support: Jensen, Dyrskjøt.
      Supervision: Dyrskjøt.
      Other: None.
      Financial disclosures: Lars Dyrskjøt certifies that all conflicts of interest, including specific financial interests and relationships and affiliations relevant to the subject matter or materials discussed in the manuscript (eg, employment/affiliation, grants or funding, consultancies, honoraria, stock ownership or options, expert testimony, royalties, or patents filed, received, or pending), are the following: Lars Dyrskjøt has sponsored research agreements with C2i, AstraZeneca, Natera, Photocure, and Ferring; has an advisory/consulting role at Ferring and UroGen; and is Chairman of the Board of BioXpedia A/S. Jørgen Bjerggaard Jensen is proctor for Intuitive Surgery; is an advisory board member for Olympus Europe, Ambu, Cepheid, and Ferring; and has sponsored research agreements with Medac, Photocure ASA, Cepheid, and Ferring. The remaining authors have nothing to disclose.
      Funding/Support and role of the sponsor: This work was funded by Independent Research Fund Denmark, The Novo Nordisk Foundation, Aarhus University (AUFF NOVA) and The Leo & Anne Albert Institute for Bladder Cancer Care and Research. The funding bodies played no direct role in the study.
      Acknowledgments: We would like to thank all the technical personnel at the Department of Molecular Medicine, Urology and Oncology, Aarhus University Hospital, for sample handling and processing. We would also like to thank GenomeDK and Aarhus University for providing computational resources and support that contributed to these research results.
      Data availability: The current legislation on data sharing requires that sensitive (including pseudonymized) data can only be shared following project approval by the National Committee on Health Research Ethics in Denmark and by the Danish Data Registry.

      Appendix A. Supplementary data

      The following are the Supplementary data to this article:

      References

        • Dyrskjøt L.
        • Thykjaer T.
        • Kruhøffer M.
        • et al.
        Identifying distinct classes of bladder carcinoma using microarrays.
        Nat Genet. 2003; 33: 90-96
        • Kamoun A.
        • de Reyniès A.
        • Allory Y.
        • et al.
        A consensus molecular classification of muscle-invasive bladder cancer.
        Eur Urol. 2020; 77: 420-433
        • Lindskrog S.V.
        • Prip F.
        • Lamy P.
        • et al.
        An integrated multi-omics analysis identifies prognostic molecular subtypes of non-muscle-invasive bladder cancer.
        Nat Commun. 2021; 12: 2301
        • Thomsen M.B.H.
        • Nordentoft I.
        • Lamy P.
        • et al.
        Comprehensive multiregional analysis of molecular heterogeneity in bladder cancer.
        Sci Rep. 2017; 7: 11702
        • Sjödahl G.
        • Eriksson P.
        • Lövgren K.
        • et al.
        Discordant molecular subtype classification in the basal-squamous subtype of bladder tumors and matched lymph-node metastases.
        Mod Pathol. 2018; 31: 1869-1881
        • Sfakianos J.P.
        • Daza J.
        • Hu Y.
        • et al.
        Epithelial plasticity can generate multi-lineage phenotypes in human and murine bladder cancers.
        Nat Commun. 2020; 11: 2540https://doi.org/10.1038/s41467-020-16162-3
        • Lai H.
        • Cheng X.
        • Liu Q.
        • et al.
        Single-cell RNA sequencing reveals the epithelial cell heterogeneity and invasive subpopulation in human bladder cancer.
        Int J Cancer. 2021; 149: 2099-2115
        • Wang H.
        • Mei Y.
        • Luo C.
        • et al.
        Single-cell analyses reveal mechanisms of cancer stem cell maintenance and epithelial-mesenchymal transition in recurrent bladder cancer.
        Clin Cancer Res. 2021; 27: 6265-6278
        • Gouin 3rd, K.H.
        • Ing N.
        • Plummer J.T.
        • et al.
        An N-Cadherin 2 expressing epithelial cell subpopulation predicts response to surgery, chemotherapy and immunotherapy in bladder cancer.
        Nat Commun. 2021; 12: 4906
        • Schmøkel S.S.
        • Nordentoft I.K.
        • Lindskrog S.V.
        • et al.
        Improved protocol for single nucleus RNA-sequencing of frozen human bladder tumor biopsies.
        Nucleus. 2023; 14: 2186686https://doi.org/10.1080/19491034.2023.2186686
        • Stuart T.
        • Butler A.
        • Hoffman P.
        • et al.
        Comprehensive integration of single-cell data.
        Cell. 2019; 177: 1888-1902.e21
      1. WISP. Weighted In Silico Pathology: a novel approach to assess intra-tumoral heterogeneity. https://cit-bioinfo.github.io/WISP/.

        • Bergenstråhle J.
        • Larsson L.
        • Lundeberg J.
        Seamless integration of image and molecular analysis for spatial transcriptomics workflows.
        BMC Genomics. 2020; 21: 482
        • Taber A.
        • Christensen E.
        • Lamy P.
        • et al.
        Molecular correlates of cisplatin-based chemotherapy response in muscle invasive bladder cancer by integrated multi-omics analysis.
        Nat Commun. 2020; 11: 4858
        • Taber A.
        • Prip F.
        • Lamy P.
        • et al.
        Immune contexture and differentiation features predict outcome in bladder cancer.
        Eur Urol Oncol. 2022; 5: 203-213
        • Strandgaard T.
        • Lindskrog S.V.
        • Nordentoft I.
        • et al.
        Elevated T-cell exhaustion and urinary tumor DNA levels are associated with bacillus Calmette-Guérin failure in patients with non-muscle-invasive bladder cancer.
        Eur Urol. 2022; 82: 646-656https://doi.org/10.1016/j.eururo.2022.09.008
        • Van de Sande B.
        • Flerin C.
        • Davie K.
        • et al.
        A scalable SCENIC workflow for single-cell gene regulatory network analysis.
        Nat Protoc. 2020; 15: 2247-2276
        • Wu S.Z.
        • Al-Eryani G.
        • Roden D.L.
        • et al.
        A single-cell and spatially resolved atlas of human breast cancers.
        Nat Genet. 2021; 53: 1334-1347
        • Joanito I.
        • Wirapati P.
        • Zhao N.
        • et al.
        Single-cell and bulk transcriptome sequencing identifies two epithelial tumor cell states and refines the consensus molecular classification of colorectal cancer.
        Nat Genet. 2022; 54: 963-975