BLCA subtypes identified based on stem cell gene sets
We collected 26 stem cell gene sets representing unique self-regenerating properties (Supplementary Table 1) and quantified the scores of 26 stem cell gene sets in each sample by ssGSEA. We used the ConsensusClusterPlus package to divide all tumor samples into k (k = 2–9) different subtypes. The CDF curve based on the consensus scores achieves the best division when k = 3. Additionally, the principal component analysis results indicated that the ssGSEA scores based on the 26 stem cell gene sets were clearly divided into 3 subtypes (Figure 2A–D). We defined these three subtypes as SCE_H, SCE_M, and SCE_L (Figure 2E). Similarly, we performed the same clustering and subtyping for the remaining data sets E-MATB-4321, GSE13507, GSE31684, GSE32548, and GSE32894.
When using ESTIMATE to evaluate the level of immune infiltration in all data sets, we found that the SCE_H immune score in all 6 data sets was much higher than that for other subtypes, and SCE_L showed the lowest immune score (Supplementary Figure 1). Comparison of the stromal content revealed the same trend (SCE_H > SCE_M > SCE_L). However, comparison of tumor purity of the BLCA stem cell subtypes showed the opposite results. SCE_H and SCE_L showed the lowest and highest tumor purity (SCE_L > SCE_M > SCE_H), respectively. This is consistent with the results observed for most HLA genes and immune cell marker genes evaluated, such as CD8A (CD8 T cells), GZMA (cytotoxic cells), IFNG (Th1 cells) PMCH (Th2 cells), CD68 (macrophages), and IL17A (Th17 cells), among others, which were significantly up-regulated and down-regulated in SCE_H and SCE_L, respectively (Supplementary Figure 2).
Based on the close relationship between BLCA stem cell subtypes and immunity, we focused on the differential expression of CD274 (PD-L1) in each subtype (Figure 3). In all six data sets, SCE_H and SCE_L showed the highest and lowest expression levels, respectively, revealing that the BLCA subtype SCE_H was more sensitive to anti-PD-1 immunotherapy than the other subtypes. Subsequent immune checkpoints inhibitor treatment response prediction and survival analysis confirmed these results.
Survival of patients with different BLCA stem cell subtypes
Because BLCA is a heterogeneous disease with a high recurrence rate and shows disease progression, exploring the relationship between subtype classification and clinical prognosis is beneficial for the prognosis assessment and corresponding clinical management of BLCA. We performed OS, RFS, and PFS analysis of the six data sets. Unexpectedly, all data sets showed consistent trends (Figure 4). SCE_H and SCE_L showed the worst and best survival in prognostic analysis (SCE_L > SCE_M > SCE_H), respectively. The P values of the log rank for the OS of TCGA, GSE31684, GSE32548, and GSE32849 were 5.631e-4, 0.038, 2.158e-4, and 8.755e-6, respectively. The P-values for the log-rank of PFS for TCGA, E-MTAB-4321, and GSE13507 were 0.004, 3.976e-9, and 7.046e-4, respectively, and the log-rank P value of TCGA RFS was 0.032.
Prediction of therapeutic response of BLCA stem cell subtypes to immune checkpoint inhibitors
Based on the above results, we further evaluated the responses of the three subtypes to immunotherapy. At present, 5 PD-1/PD-L1 immunotherapy drugs have been approved by the Food and Drug Administration for treating BLCA, including nivolumab and pembrolizumab (both PD-1 inhibitors) approved in 2016 and 2017 for treating patients with locally advanced or metastatic urothelial cancer who were administered first-line platinum-containing chemotherapy for one year [24-29]. We used the tumor immune dysfunction and exclusion algorithm to predict the likelihood of a response to immunotherapy. The results showed obvious differences in the responses to immunotherapy among the SCE_H (20%, 32/158), SCE_M (42%, 71/168), and SCE_L groups (61%, 54/88) (P = 2.951e-10). We further performed subclass mapping to compare the expression profiles of the three stem cell subtypes which were defined using another published data set containing 47 patients with melanoma who responded to immunotherapy [30]. In pairwise comparison of the three subtypes, more promising results were observed in SCE_H for the anti-PD1 and anti-CTLA4 treatments compared to the other subtypes (Figure 5A–C) (anti-PD1 therapy: SCE_H vs SCE_L, FDR = 0.036; SCE_H vs SCE_M, p = 0.046; SCE_M vs SCE_L, FDR = 0.048; anti-CTLA4 therapy: SCE_H vs SCE_L, FDR = 0.036; SCE_H vs SCE_M, FDR = 0.008). We further correlated the BLCA stem cell typing results with the published molecular typing and immunotyping results in TCGA cohort. SCE_H mainly corresponded to the molecular subtypes luminal-infiltrated and basal squamous and C1 and C2 for immune subtypes; SCE_L mainly corresponded to luminal-papillary and C1–C4, and SCE_M showed a wide distribution of molecular and immune subtypes (Figure 5D).
Differences among 22 human immune cell subgroups of BLCA stem cell subtypes in CIBERSORT
To explore the reasons for the difference in survival of patients with different BLCA stem cell subtypes, we used the CIBERSORT algorithm to calculate the proportions of 22 immune cells in each subtype of the six data sets, with P <0 .05 as the threshold for screening. The results showed that the proportions of macrophages M0, M1, and M2 showed an upward trend in the SCE_H subtype (except for GSE13507) and the proportion of regulatory T cells (Tregs) was significantly (P <0 .05) increased in the SCE_H subtype (Figure 6A–F). We also used the ssGSEA scores of immune cells in each cohort as continuous variables and performed univariate Cox analysis (Supplementary Table 2). Further, we divided the median value of the ssGSEA scores of the corresponding data set into groups with high and low scores. A high score indicated that patients with a high macrophage M0 content had a worse prognosis. This is consistent with the poor clinical prognosis of patients with the SCE_H subtype compared to that of patients with the other subtypes. The trend for macrophages M2 was similar to that of macrophages M0, whereas macrophages M2 and Tregs showed opposite trends (Figure 6G–L). This reveals that compared to Tregs, macrophages M0 and M2 have completely opposite regulatory mechanisms in the prognosis of BLCA.
Correlation of CD274 with stemness genes and risk observation of stem cell subtype populations
We identified an important immune role for CD274 in the stem cell subtypes, and thus we further explored the correlation between CD274 and the identified stemness genes. Figures 7A–G show scatter plots of the expression of CD274 and stemness genes CD44, GATA3, HIF1A, ID1, MYC, SOX9, and CXCL8 in BLCA TCGA cohort. Among them, CD274 was negatively correlated with ID1 and GATA3 and positively correlated with CD44, HIF1A, MYC, SOX9, and CXCL8 (Figure 7H). We divided the patients according to the optimal expression cut-off of CD274 and each stemness gene into risk groups I, II, III, and IV (for example, CD274 and CD44 corresponded to CD274lowCD44low, CD274highCD44low, CD274low CD44high, and CD274highCD44high). According to the scatter plot, for each pair of risk groups divided by CD274 and the optimal threshold of stemness gene expression, patients with stem cell subtypes were mainly concentrated in groups I and III (87–88%), the SCE_L subtype was mainly concentrated in risk group III of the CD274 and GATA3 and ID1 pairs, and in the CD274 and CD44, HIF1A, MYC, SOX, and CXCL8 gene pairs were concentrated in risk group I. Further survival analysis of groups I and III of each gene pair showed that patients with the higher SCE_L subtype population had better survival than patients with a lower SCE_L subtype population (Figure 7I–O). This is consistent with observation that patients in the SCE_L group had the longest survival time.
Stemness index based on mRNA expression (mRNAsi) in BLCA
We evaluated the differences in stem cell subtypes defined based on the stemness index obtained from an existent one-class logistic regression (OCLR). Five of the six data sets (Figure 8A–E) showed the decreasing trend of SCE_M > SCE_L > SCE_H, which showed that different levels of stem cells enriched the potential for self-renewal and differentiation. In the E-MTAB-4321 data set with <T2 accounting for 96.6% (460 cases), the opposite results were observed. The stemness index decreased according to the trend of SCE_H > SCE_L > SCE_M (Figure 8F), reflecting the high impact of the T stage on the stemness index and revealing dynamic development of the tumor. This difference was also reflected in normal and cancer tissues, primary and recurrent (Figure 8G–I); thus, the stemness index of cancer tissue is much higher than that of normal tissue, with recurring cancer showing a higher index than primary cancer. Even in the non-uniform molecular classification of BLCA, the stemness index can effectively characterize the characteristics of different subtypes. According to the BLCA subtypes described by Robertson, neuronal and luminal infiltrated subtypes have the highest and lowest stemness indices respectively, which is consistent with corresponding reports showing that the neuronal type has the worst prognosis compared to other types and a better prognosis for the luminal infiltrated type. Similarly, among the five BLCA subtypes classified by Lund University [8], the genomically unstable and infiltrated types have the highest and lowest stemness indices, respectively (Figure 8J–K). Figures 8L–P show the difference in the T stage stemness index in each data set.
GSEA for BLCA stem cell subtypes
To explore the biological changes caused by differences in the enrichment of stem cells, we conducted pairwise comparison of the GSEA results for each subtype.
By selecting at least one pathway with an FDR <0.05 for the same pathway, we found that as the enrichment of stem cells increased, epithelial mesenchymal transition (EMT) became more significant. EMT is considered as a signal of malignant transformation in all cancers, giving cells the ability to metastasize and invade, including by imparting stem cell characteristics, reducing apoptosis and aging, and resisting chemical and immunotherapy [31, 32] (Table 1). EMT can also activate multiple pathways, regulate cell metabolism, angiogenesis, proliferation, and migration, and enable cells to respond hypoxic environments. Pathways are also significantly enriched, for example, in hypoxia, angiogenesis, inflammatory response, IL6-mediated Jak-STAT signaling pathway, and KRAS signal up-regulation (Figure 9). These pathways together constitute a vicious circle of cancer occurrence, proliferation, invasion, and metastasis.
Table 1
GSEA for BLCA stem cell subtypes.
pathway | NES | FDR |
HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION | | |
SCE_H vs SCE_L | 2.21 | 0.005 |
SCE_H vs SCE_M | 2.42 | 0 |
SCE_M vs SCE_L | 2 | 0.008 |
HALLMARK_HYPOXIA | | |
SCE_H vs SCE_L | 1.8 | 0.031 |
SCE_H vs SCE_M | 1.41 | 0.129 |
SCE_M vs SCE_L | 1.72 | 0.034 |
HALLMARK_INFLAMMATORY_RESPONSE | | |
SCE_H vs SCE_L | 2.01 | 0.018 |
SCE_H vs SCE_M | 1.81 | 0.019 |
SCE_M vs SCE_L | 1.89 | 0.015 |
HALLMARK_IL6_JAK_STAT3_SIGNALING | | |
SCE_H vs SCE_L | 1.84 | 0.034 |
SCE_H vs SCE_M | 1.61 | 0.069 |
SCE_M vs SCE_L | 1.71 | 0.033 |
HALLMARK_ANGIOGENESIS | | |
SCE_H vs SCE_L | 1.76 | 0.037 |
SCE_H vs SCE_M | 1.98 | 0.004 |
SCE_M vs SCE_L | 1.33 | 0.182 |
HALLMARK_KRAS_SIGNALING_UP | | |
SCE_H vs SCE_L | 1.94 | 0.019 |
SCE_H vs SCE_M | 2.19 | 0.001 |
SCE_M vs SCE_L | 1.50 | 0.087 |
NES: normalized enrichment score; FDR: false discovery rate. NES > 1.0 and FDR < 0.25 was considered to be significant. |
Key gene networks identified in BLCA stem cell subtypes
A total of 14 gene modules were generated based on WGCNA. Among them, brown and black modules showed the strongest relationships with stem cell subtypes, and the brown module was positively correlated with SCE_L and negatively correlated with SCE_H and SCE_M. The black module was positively correlated with SCE_H and negatively correlated with SCE_M and SCE_L (Figure 10A). Two pairs of transcription factor regulators were identified in the brown and black modules (Figure 10B–C), IRF5 and GATA3 and GRHL2 and GATA6. IRF5 is a key transcription factor regulating the differentiation of M1 macrophages into M2 macrophages, enabling M2 macrophages to play an anti-inflammatory role; these cells also have an important role in tissue repair and reconstruction and cancer occurrence [33-35]. GATA3 is type 2 helper T cell (Th2) cytokine-specific transcription factor and a key stemness gene that regulates cell differentiation. It enables Th2 to express IL-4 and other cytokines, promotes antibody production, mediates humoral immunity, and suppresses anti-tumor immunity [36-38]. Therefore, GATA3 may act as a tumor suppressor gene in BLCA. GRHL2 and GATA6 play important regulatory roles in many life activities such as embryonic development, damage repair, epidermal barrier formation, tracheal epithelial formation, and neural tube development, but they also play an important role in the occurrence and development of tumors, cell proliferation, invasion, and metastasis are closely related, which were verified in the relevant pathways identified in the black module [39-45] (Figure 10D–E). Thus, we identified GATA3 and GATA6 as important transcription factors with opposite expression and effects on prognosis in patients with different BLCA stem cell subtypes (Figure 10F–I).