Development and validation of an intra-tumor heterogeneity-related signature to predict prognosis of bladder cancer: a study based on single-cell RNA-seq

Abstract

Intra-tumor heterogeneity (ITH) was a potential mechanism of progression and drug resistance in bladder cancer (BCa). However, the understanding of ITH in BCa remains insufficient. Single-cell RNA sequencing (scRNA-seq) profiles of 2075 cells were analyzed, and 2940 cell markers were screened. The ITH of 396 cases was evaluated, and 96 ITH-related genes were identified. Based on the gene-pair strategy, 96 genes were cyclically paired, and an 8-gene-pair model was successfully established to evaluate the overall survival of BCa through Lasso and multivariate Cox regressions. The risk model showed high predictive value in the training dataset (n = 396, p = 0) and external validation datasets (n = 165, p = 2.497e-02; n = 224, p = 3.423e-02). The model was also valuable for the prediction of clinical treatment outcomes. Totally, a prognostic model based on scRNA-seq and ITH was successfully constructed and validated in large cohorts, providing novel clues for ITH study of BCa.

Keywords: bladder cancer; intra-tumor heterogeneity; prognosis; risk model; scRNA-seq.

Figure 1

Evaluation of ITH with MATH in BCa. (A) The workflow of this study. (B) The patients with high MATH values suffered an unfavorable prognosis. (C) The estimated MATH values acted as a potential predictor for chemosensitivity with Wilcoxon signer-rank test. (D) MATH values were positively correlated with TMB. (E) The cases in the high-MATH group had significantly lower immune and stromal components in TME. (F) Spearman correlation analysis indicated MATH values were negatively correlated with routine immune checkpoint genes, including PD1, LAG3, GAL9, CTLA4, TIM-3, and TIGIT. The red lines and green lines represented positive correlation and negative correlation, respectively. The boldness of the lines was positively associated with the strength of the correlation. ^*p < 0.05; ^**p < 0.01; ^***p < 0.001.

Figure 2

Characterization of scRNA-seq from 2075 cells. (A) Quality control plots of cell samples. (B) The sequencing depth was negatively correlated with the proportion of mitochondrial genes and positively associated with detected gene numbers. (C) 1500 variable genes across cell samples were identified. (D, E) PCA was conducted to reduce the dimension of data sets. (F) Cell samples were classified into 14 clusters with the t-SNE algorithm. (G) The trajectory analysis of 14 cell clusters.

Figure 3

Establishment of HRS. (A) 96 of 2075 cell marker genes were significantly correlated with MATH. (B, C) GO functional annotation (B) and KEGG pathway enrichment (C) of the 96 genes. (D, E) 10 crucial gene-pairs correlated with OS were identified via Lasso-Cox regression. (F) 8 of 10 gene pairs were included in the prognostic model by multivariate Cox analysis with stepwise. (G) The expression values of 13 genes comprising the 8-gene-pair signature in 14 cell subpopulations. (H) The mutational rates of the 13 genes in different stages in BCa. The size of the bubble represented the mutational rates in all samples.

Figure 4

Validation of HRS. (A) A nomogram was plotted to visualize the HRS. (B, C) The calibration curves for 3- (B) and 5-year (C) OS prediction. (D–F) Kaplan-Meier survival analysis of HRS in TCGA (D), GSE13507 (E) and GSE32894 (F) cohorts. (G–I) The distribution of HRS in TCGA (G), GSE13507 (H) and GSE32894 (I). (J–L) The distribution of survival status in TCGA (J), GSE13507 (K) and GSE32894 (L).

Figure 5

HRS showed superiority in OS prediction compared with the clinicopathological features. (A, B) HRS was an independent prognosis predictor in univariate (A) and multivariate (B) analyses. (C–E) ROC analysis indicated HRS had better ability than the clinical traits in 1- (C), 3- (D) and 5-year (E) OS prediction.

Figure 6

The correlation between HRS and other clinicopathological variables. (A) The strip curve displayed HRS was significantly correlated with pathological tumor stages. (B–G) Correlation analysis of HRS with age (B), gender (C), pathological stages (D), pathological T stages (E), pathological N stages (F), and M stages (G).

Figure 7

HRS was a potential predictor for clinical treatment of BCa. (A, B) HRS was significantly correlated with MATH, codetermined by difference (A) and correlation (B) analysis. (C) GAL-9 was differentially expressed between low- and high-HRS groups. (D) The cases in the high-HRS group were more likely to be associated with the high infiltration of M0 macrophages, activated mast cells, and neutrophils, whereas they were negatively correlated with the infiltration of CD8⁺ T cells and Tregs. (E) High HRS was linked to a lower IC50 for chemotherapeutics like cisplatin and doxorubicin, whereas it was correlated to a higher IC50 of methotrexate. ^*p < 0.05; ^**p < 0.01; ^***p < 0.001.

Figure 8

Functional enrichment analysis. (A) The heatmap showing the analysis results of GSVA. (B, C) Four pathways, including G2M checkpoint, mitotic spindle, mTORC1 signaling, and epithelial mesenchymal transition, were positively related to HRS, which was codetermined by GSEA and GSVA. (D, E) Two pathways, including DNA repair and oxidative phosphorylation, were negatively associated with HRS through the combined analysis by GSEA and GSVA.

Figure 9

The candidate compounds targeting HRS. (A) The volcano plot displayed the DEGs between low-HRS and high-HRS cases. (B–I) Seven compounds, including cephaeline (C), LY-294002 (D), lycorine (E), naltrexone (F), nefopam (G), tanespimycin (H), and wortmannin (I), were identified.