Introduction Of Main Rcd Types And The Workflow In This Study
As shown in Fig. 1A, the main morphological characterization of apoptosis is cell shrinkage, nuclear condensation, karyorrhexis, and apoptotic body formation23. Unlike apoptosis, pyroptosis could be triggered by activation of CASP1 or CASP11, the Gasdermin D (GSDMD) formed a GSDMD-C and GSDMD-N fragment. The GSDMD-N binds to cell membrane phospholipids and induces the formation of cleavage pores, which leads to membrane lysis23–26. Ferroptosis is driven by iron accumulation-induced ROS production and lipid peroxidation, which lead to dysmorphic small mitochondria with decreased crista, and condensed, ruptured outer membranes23. In the process of autophagy, the main steps include the fusion of the autophagosome with the lysosome and autolysosome, and membrane rupture27. Necroptosis is mainly regulated by activation of threonine kinase 3 (RIPK3), activation of mixed lineage kinase domain-like pseudokinase (MLKL), and membrane rupture28. These RCD mentioned above were the primary manner of programmed cell death, most of which are related to the rupture of the plasma membrane10. To explore the cross-talk, molecular mechanisms, transcriptional regulation, and clinical relevance of RCDs in bladder cancer, we generated a 5 step analysis pipeline: i, Checking databases and relevant literature to determine these five types of RCD-related genes, then reclassified these genes, identifying the complex relationship is exist in mRNA correlation and PPI network analysis; ii, Analyzing the cross-talk and regulatory relationship of RCDGs in the aspect of somatic mutation, CNVs, DNA methylation, RNA modification, ceRNA network, and regulator network; iii, based on the gene Dependency Score (gDS), IC50 of chemotherapeutic drug, aneuploidy score, and radiosensitivity index, investigating the importance of RCDGs in bladder cancer; iv, RCD cluster construction, molecular mechanisms, immune infiltration, and therapy; v, exploring the relationship between RCD clusters and consensus cluster, and the distribution of RCDGs in BLCA single-cell subtypes (Fig. 1B). Referring to the public database and high-quality literature, we ultimately identified 445 apoptosis genes, 254 ferroptosis genes, 222 autophagy genes, 40 pyroptosis genes, and 28 necroptosis genes (Fig. 1C, Table S1a). Venn diagram showed that several genes functioned in different RCD types (Fig. 1D, Table 1). To facilitate subsequent studies, we named these genes as multifunctional genes because they do not fall into any single category of RCD.
Cross-talk Of Rcdgs In Bladder Cancer
To determine the underlying cross-talk of RCD, we calculated the Pearson Correlation Coefficient (PCC) based on the expression of RCDGs (Fig. S1A, details in Table S2a). Unsupervised clustering was performed for the correlation coefficients between RCDGs’ expression, and we divided the correlation coefficients diagram into five zones (Fig. S1A, Table S2b). We found that the genes in zone1, zone3, zone3, and zone5 showed a significant positive or negative correlation. The genes in zone4 showed a lower correlation with each other or other zones’ genes. Interestingly, genes in different regions functioned in different types of RCD, not a specific type of RCD gene is distributed in a certain region (Fig. S1B), revealing that co-expression exists in different RCDGs. We further analyzed the protein-protein interaction of these RCDGs and visualized them in a network diagram. We found an extremely complex PPI network, and most genes representing different types of RCD are involved in this network (Fig. S1C, Table S3). These results revealed that RCDGs have complex relationships.
Mutation Characterization Of Rcdgs
To explore the mutation characterization of RCDGs in bladder cancer, we downloaded the mutation data from the public database, including TCGA-BLCA, ICGC-BLCA, and CCLE-BLCA datasets. In the TCGA-BLCA dataset, we found that the most common variant type of RCDGs is missense mutations, SNP, and C > T alterations (Fig. 2A). We randomly selected 1000 genes from the TCGA-BLCA mutation data and calculated the mutation rate of random genes and RCDGs. We found mutation rate of the genes representing different RCD types was higher than that of random genes, especially pyroptosis and multifunction genes, indicating that RCDGs occurred with larger mutation disturbance (Fig. 2B, details in Table S4). Among these, 86 RCDGs showed a higher mutation rate (> 2.5%) than other RCDGs in TCGA-BLCA (Fig. 2C). The mutation rate of TP53, PIK3CA, RB1, ATM, BIRC6, CREBBP, ERBB2, SPTAN1, and ERBB3 was over 10%. Some of these gene mutations influence patients' prognosis and immune therapy response29–31. What is more, several RCDGs’ expressions were affected by somatic mutation. In wild type, the expression of BAP1, PTEN, RB1, SPTAN1, TP53, and TSC1 were lower than that of mutation type. Whereas CDKN2A, ERBB3, HRAS, NFE2L2, and NLRP3, their expression was higher than that of mutation type (Fig. S2A). We further calculate the mutation rate using CCLE-BLCA and ICGC-BLCA datasets (Fig. S2B and C). TP53, PIK3CA, CREBBP, RB1, and ERBB3 also showed high mutation frequency in these two datasets. Interestingly, among the top 25 high mutation rate genes of the TCGA-BLCA dataset, we found that seven genes were multifunction RCDGs (P-value = 0.016; Hypergeometric test) (Fig. 2D). Among the top 25 high mutation rate genes of the ICGC-BLCA dataset, seven genes were also multifunction RCDGs (P-value = 0.016; Hypergeometric test) (Fig. S2D). Between the top 25 mutation genes of TCGA-/ICGC-BLCA dataset, 10 genes showed common higher mutation rate (P-value = 0.047; Hypergeometric test) (Fig. S2E), 5 of them were multifunction RCDGs (P-value = 0.0068; Hypergeometric test) (Fig. 2E). These results reveal that multifunction genes contribute more to the tumor mutation load than other genes. Other RCDGs with a high mutation frequency may play roles in different cell death, which need further validation. For ICGC and TCGA, five multifunction genes showed the highest mutation, including TP53, RB1, PIK3CA, ERBB2, and ATM (Fig. 2F); they also interacted tightly in protein levels. E545K/Q/A, R248Q/W/G, X405_splice, E2164Q/G, and S310F/Y were the hotspot mutation of PIK3CA, TP53, RB1, ATM, and ERBB2 in BLCA, respectively (Fig. S2F). Some of these hotspot mutations have been reported to influence BLCA or other cancer types32–34. We further analyzed the cross-talk of the 86 high mutation rate genes of TCGA-BLCA. Most RCDGs showed mutation co-occurrence (Fig. 2G and Fig. S2G, details in Table S5). For example, TP53 (Multifunction genes) mutation may lead to the mutation of NLRP9 (Pyroptosis gene), CDKN2A (Multifunction genes), and RB1 (Multifunction genes); DPYD (Apoptosis genes) mutation may co-occur with CEBPZ (Apoptosis genes), ABCC1 (Ferroptosis genes), TRPM7 (Necroptosis genes), PLCG1 (Pyroptosis genes). In contrast, PIK3CA (Multifunction genes) with PSME4 (Apoptosis genes), RB1 (Multifunction genes) with NFE2L2 (Multifunction genes) showed mutually exclusive mutation.
Copy Number Characterization Of Rcdgs
Copy number variants (CNVs) are also an essential source of genetic diversity and can drive rapid adaptive evolution and progression of heritable and somatic human diseases, such as cancer35. To explore the CNV characterization of RCDGs in bladder cancer, we calculated the CNV rate of each RCDGs in TCGA-BLCA. We found that all RCDGs of chromosome 8 and several RCDGs of chromosomes 9, 17, and 20 showed a higher CNV rate (> 60%) (Fig. 3A and B, Table S6a). Similar results were observed in CCLE-BLCA (Fig. S3A, Table S6b). Among these, 43 RCDGs showed a higher CNV rate in both datasets (P-value = 0, Hypergeometric test) (Fig. S3B). Furthermore, we found that the pyroptosis genes showed a higher copy number amplification rate, whereas multifunction genes showed a higher copy number deletion rate (Fig. 3C, Fig. S3C, and D). For the high mutation rate genes of TCGA-BLCA, several genes (including PIK3CA, ETBB2, ERBB3, and BIRC6) occurred to a higher copy number amplification rate, whereas TP53, ATM, RB1, and CERBBP, showed a higher copy number deletion rate (Fig. 3D, Table S6c). We further depicted the landscape of RCDGs with amplification or deletion rates > 50% (Fig. 3E and F; Fig. S3E and F). 43 RCDGs showed a consistent high amplification rate in these two datasets, and these genes distribute into different RCD types (P-value = 0, Hypergeometric test) (Fig. 3G). We also observed 23 genes with a consistently high deletion rate in two datasets (P-value = 0, Hypergeometric test) (Fig. 3H). Interestingly, we observed that the CNV types of several genes showed highly consistent in BLCA, such as BLC2L1, CHMP4B, and E2F1; TFRC and PAK2, etc., which seems to imply that copy-number variants of these genes may co-occur in BLCA (Fig. 3E and F; Fig. S3E and F). Therefore, we further calculated the correlation between CNVs type of these high amplification or deletion rate genes. Multiple genes showed a high correlation with each other (Fig. 5I and J). Notably, such as SRC and BLCAP, PKT6 and RGS19, TFRC and PAK2, RPL8 and HSF1, showed a completely consistent CNV charaterization in the TCGA-BLCA (Spearman r = 1). These genes play roles in different RCD types and showed high co-occurrence probability in DNA copy number variate levels. Copy numbers not only drive human disease but also can influence gene expression for the majority of genes36. In our study, we also found that multiple types of RCDGs were highly correlated (PCC > 0.3, adj.Pvalue < 0.05) with their linear copy number levels (Fig. S3G, Table S6d). The expression disturbance of these genes may be the product of their copy number variants.
Dna Methylation Regulation Of Rcdgs
DNA methylation is an epigenetic mechanism closely correlated to cancer occurrence and development37. DNA methylation regulates gene expression by recruiting proteins involved in gene repression or inhibiting transcription factors' binding to DNA38. To explore RCDGs driven by DNA methylation, we utilized the ChAMP R-package to analyze the DNA methylation characterization of RCDGs in BLCA. The site of DNA methylation of RCDGs is widely distributed in chromosomes 1–22 and X, especially in chromosomes 1 and 17 (Fig. 4A). For different genomic regions, the site of DNA methylation is mainly located in the opensea and island (Fig. 4B). From the point of view of gene distribution, DNA methylation sites of RCDGs were mainly distributed in the body and promoter region (including TSS200, TSS1500, 5'UTR, and 1stRxon) (Fig. 4C). We further calculated the differential methylated sites (DMSs) of RCDGs comparing bladder cancer with normal control (|LogFC|>0.2 & adj.P.value < 0.05). There were 181 hypermethylated and 493 Hypomethylated DMSs in BRCA (Fig. 4D, Table S7a), which correspond to 72 hypermethylation and 218 hypomethylated genes. The methylation levels of these DMGs could district the tumor and normal tissues (Fig. 4E). Among these, the hypermethylation sites were mainly located on the island, and the hypomethylation sites were mainly in opensea (Fig. 4F). From the view of gene distribution, the hypermethylation sites were mainly located in TSS200 and the body, whereas the hypomethylation sites were mainly located in Body TSS1500 (Fig. 4G). To further understand which genes are driven by DNA methylation, the Pearson correlation test was performed on these DMSs to find the DMSs that can drive gene expression (|PCC|> 0.3 & adj. P-value < 0.05) (Table S7b). 109 DMGs sites were negatively correlated with their correspondent genes, and 151 DMGs sites were positively correlated with their correspondent genes. Among these, the expression of TP63, IRAK3, ERBB3, and BCAT2 showed a highly negative correlation with DNA methylation levels (|PCC| > 0.6 & adj.Pvlaue < 0.05), revealing that they were driven by DNA methylation. Interestingly, for the negative correlation relationship, the main genes were apoptosis and ferroptosis genes, and the main location was the body and promoter (including TSS200, TSS1500, 5'UTR, and 1stExon) region (Fig. 4H and I). In contrast, the main positive genes were apoptosis and autophagy genes, and the main location was the body and the 3'UTR region (Fig. 4J and K).
Rna Modification Of Rcdgs
Accumulate evidence indicated that RNA modification pathways are also misregulated in human cancers and may be ideal targets for cancer therapy39. The main RNA modification is m6A and m5C modification, which were widespread epigenetic modifications in tumorigenesis40–42. The RNA modification is determined by the coordinated actions of the three regulators: methyltransferases (writers), RNA-binding proteins (readers), and demethylases (erasers)43. The main m6A regulators which affect mRNA are YTHDF1, YTHDF2, YTHDF3, YTHDC1, YTHDC2, METTL3, FTO, and ALKBH5 (Table 2). Through the Pearson correlation analysis, we observed that multiple RCDGs were highly positively or negatively correlated with these m6A regulators (Fig. S4, Table S8). Interestingly, few protein-protein interactions were discovered between m6A regulators and RCDGs, indicating that the regulation region of these m6A regulators on RCDGs was mainly not in the protein region but most likely in the mRNA region. However, when we wanted to study further the influence of m6A level on mRNA level of RCDGs in bladder cancer, there was no data for utilizing in public dataset. A recent study involving the m5C methylation profiling by high throughput sequencing in BLCA provided a foundation for exploring the RCDGs which might be regulated by m5C methylation44. Chen et al. identified 5221 m5C sites with differential significant methylation levels, and this site represents 1902 genes (Table S9a). Among these, 266 m5C sites representing 84 RCDGs showed differential methylation expression, including 56 hypomethylation sites and 210 hypermethylation sites (Fig. 5A and B, Table S9b). These genes mainly come from apoptosis, autophagy, ferroptosis, multifunction, and pyroptosis (Fig. 5C). The m5C levels of the genes influenced the apoptosis (DCN, NGF, RHOT2, DIABLO, and LMNA), autophagy (AMBRA1, RPTOR, and UVRAG), Ferroptosis (TXNRD1, AKR1C3, AKR1C2, and YWHAE), and multifunction (GZMB) were downregulated in bladder cancer. In contrast, the m5C levels of other RCDGs were upregulated in bladder cancer (Fig. 5C). In addition, Chen et al. also detected the expression profiling by high throughput sequencing in BLCA and adjacent normal tissues. Through a differential expression analysis, we found that the mRNA expression of several RCDGs was upregulated or downregulated (Fig. 5D, Table S9c). We further calculated the correlation between these DEGs and their correspondent m5C levels, and the results showed that the mRNA expression of several genes might be regulated by their m5C levels (Table S9d). For example, TP73 was upregulated in mRNA, and m5C methylation levels in BLCA, the mRNA expression and m5C levels of it show a positive correlation (PCC = 0.3553, P-value = 0.0335) (Fig. 5E); DCN also downregulated in mRNA and m5C levels in BLCA, showed a positive correlation between mRNA expression and m5C levels (Table S9d); PIK3R3 and GZMB, were downregulated in mRNA levels and upregulated in m5C levels, showed a negative correlation between mRNA expression and m5C levels, etc. Combining previous reports that m5C levels could influence the mRNA stability44, it is considered that these mRNA expressions of RCDGs were affected by their m5C levels. Of note, there is heterogeneity among samples and different sequencing platforms, and the results often differ. The m5C regulator may regulate several genes, but the alteration of m5C levels has not reached the thresholds of DMGs. Therefore, we further analyzed the relationship between the main mRNA-related m5C regulators regulated and RCDGs (Fig. 5F). We found that multiple RCDGs were positively or negatively correlated with these regulators (|PCC| > 0.3, adj. Pvalue < 0.05) (Fig. S5, Table S8). Like m6A regulators, most RCDGs did not interact with m5C regulators in the protein layer. It is considered that the regulators may regulate these genes in an m5C-dependent manner. A recent study showed that ALYREF could regulate PKM2 in an m5C-dependent manner. In our results, the correlation between mRNA expression of ALYREF and PKM2 was lower in two public BLCA data (Fig. 5G and H). However, the ALYREF RIP-seq data showed that m5C sites of PKM2 could interact with ALYREF45. Therefore, we also tried to find the m5C site of RCDGs which interacted with ALYREF. From the Yang et al. ALYREF RIP-seq data, we found that multiple types of RCDGs could interact with ALYREF in their m5C site, and some of them showed a positive correlation with ALYREF in mRNA levels (Fig. S6A, Table S9e), such as ATG3, ATG13, ATM, etc. These results further provide evidence for us to identify RCDGs regulated by ALYREF in m5C dependent manner. We further analyzed a high throughput sequencing data - GSE133621, including the T24-siNC and T24-siNSUN2 RNA-seq data. Through a differential expression analysis, we identified 223 RCD DEGs, including 158 upregulated genes and 65 downregulated genes (Fig. S6B and C, Table S10). Several down-regulated genes showed a positive correlation with NSUN2 in mRNA levels, such as CDK2, PSMD3, RAB7A, USP10, etc., and they do not interact with NSUN2 in protein levels, revealing their mRNA transcription may be regulated by NSUN2. In GSE175714 data, we analyzed the DEGs in sgYBX1, and sgLuci data, including 52 RCD upregulated RCDGs, and 73 downregulated RCDGs (Fig. S6D and E, Table S11). Several down-regulated genes positively correlated with YBX1 in mRNA levels, such as PSMF1, GAPDH, and MYC, revealing that their mRNA transcription may be regulated by YBX1. However, whether these regulatory mechanisms are dependent on m5C requires further experimental verification.
Transcriptional Regulation Network Of Rcdgs
Control of gene and protein expression is required for cellular homeostasis and is disrupted in disease46. Micro-RNA are small non-coding RNAs that play a wide variety of critical roles in post-transcriptionally regulating gene expression47. Competitive endogenous RNAs (ceRNAs) could sequester miRNA and diminish its repressive effects on other transcripts46. To explore the ceRNA regulatory network of RCDGs, we firstly analyzed the expression of RCDGs and lncRNA in BLCA. The interactions with |PCC| > 0.6 and adjust. P-value < 0.05 are shown in Fig. 6A (Table S12a). Meanwhile, we downloaded the highly conserved microRNA families from the mircode database to predict the potential miRNA competition with lncRNA (Table S12b). A ceRNA network was established, and the results showed a complex regulatory relationship between lncRNA-miRNA-RCDGs (Fig. 6A, Table S12c). In addition, previous research reported that 23 regulator genes were associated with bladder cancer48, so we further assessed the mRNA correlation between 23 regulator genes and RCDGs, the interaction with |PCC| > 0.6, and adjusted. P-value < 0.05 are shown in Fig. 6B (Table S13). RCDGs are highly associated with multiple regulators, revealing that these regulators may regulate RCDGs in RNA expression.
The Importance Or Therapeutic Relevance Of Rcdgs
To understand the importance of these RCDGs, we firstly utilized a gene dependency score (gDS) to estimate the “Essentiality” of RCDGs in bladder cancer cells49 (Table S14a). The genes with gDS < -1 indicated their deletion severely affects cell viability; gDS between − 1 and 0 indicated these genes influence cell viability to some degree but are not lethal; gDS > 0 of genes indicated whose deletion does not affect cell viability. Among these, 214 genes with gDS < 0 in at least one of the three bladder cancer cells, we called these genes potential essential genes (Fig. 7A, Table S14b). Ninety-three genes with gDS > 0 in all the three bladder cancer cells, we called these no essential RCDGs (Fig. S7A, Table S14c). Among these, the deletion of several RCD-promoted genes showed no essential for cell viability, such as CASP2, RIPK3, TFRC, GLS2, and ACSL4, indicating that up-regulated these genes could kill tumor cells in different RCD manner. In contrast, the deletion of several RCD-inhibited genes showed essential for cell viability, such as NFE2L2, AURKA, BCL2L2, and BCL2L10, indicating that up-regulated these genes could play a cancer-promotion role by inhibiting their correspondent RCD types. Aneuploidy, indicating the copy number alteration of whole chromosome arms or chromosomes, has been revealed to be a general and remarkable feature of malignancy genomes50. For the RCDGs, we found several genes (including TOP2A, MSH2, and E2F1) showed a positive correlation with the aneuploidy score, whereas NLRP1 and GNB2L1 showed a negative correlation with the aneuploidy score, indicating these genes could drive the aneuploidy in BLCA (Fig. S7B, Table S15).
For advanced bladder cancer, neo-adjuvant therapy and traditional chemotherapy were mainly therapeutic options. Current chemotherapy therapy drugs including gemcitabine51, cisplatin51, methotrexate52, vinblastine52, paclitaxel51. Drug sensitivity influenced the therapeutic effect. To determine drug sensitivity-related genes, we calculated the estimated IC50 of main drugs in CCLE BLCA cells and TCGA- BLCA tissues; the results were consistent in these two datasets (Fig. S7C and D, Table S16a and b). Then we calculated the PCC between RCDGs’ expression and IC50 of the drugs, and the correlation network with |PCC| > 0.3 and adj.P.value < 0.05 were identified (Fig. 7B and Fig. S7E, Table S16c and d). The positive correlation indicated that the genes might increase the IC50 of the drug, and the negative correlation indicated that the genes might decrease the IC50 of the drug. Though the CCLE samples are small, multiple correlations did not reach the threshold, but we also observed some consistent results. For example, TOP2A could decrease the IC50 of the drugs, revealing it may play a potentially important role in chemotherapeutic sensitization.
In addition to cystectomy and chemotherapy, radiotherapy (RT) for bladder cancer is an effective alternative53. To understand the effect of RCDGs’ alteration for radiotherapy, we used a radiosensitivity index (RSI) constructed by a previous study54. High RSI indicated radioresistance, and Low RSI indicated no radioresistance54,55. In the TCGA-BLCA cohort, we found that the prognosis of the patients in the high RSI group is poor than that in the low RSI group (Fig. S7F). RSI was an independent prognostic factor for BLCA patients (Fig. S7G). In addition, in the cohort of patients receiving radiation therapy, the prognostic of the high RSI group was worse (Fig. 7C). These results revealed that RSI could predict the prognosis and radiotherapy effect of BLCA patients. Further analysis revealed that several RCDGs (including PSME1, IKBKB, ALOX5, BAG1, GSDMB, PSMB10, IRF1, and IKBKE) were negatively correlated with RSI, representing that these genes could enhance the radial sensitivity (Fig. 7D, Table S17a and b).
Clinical Relevance And Molecular Mechanism Of Rcd Clusters
To further explore the effect of complex RCD relationships on clinical outcomes, we first conducted the TSNE analysis (Fig. 8A). On a two-dimensional axis, the normal sample straddles the tumor sample, indicating that the expression of the RCD gene alters in two different directions. Interestingly, the BLCA patients were divided into two clusters based on the expression of RCDGs, and the prognosis of the patients in cluster1/cluster2 was significantly discriminative (Fig. 8B and 8C). We conducted the same analysis in the other three MIBC datasets (including GSE13507, GSE48075, and GSE48276). The results also showed that the patients could be separated into two different OS rate groups based on the RCDGs (Fig. S8A-I). Intriguingly, when we divided the patients based on the differential expression RCDGs and prognosis-related RCDGs, the patients in clusters showed no significant difference in overall survival, revealing that all RCDGs caused the patients' prognosis of different clusters. In other cancer types, we conducted the same analysis (Fig. S9 and Fig. S10). RCD clusters of multiple cancer types showed different prognoses, revealing that the various expression characterization of RCDGs in cancers could cause different clinical outcomes. In these two clusters, we observed that the number of C > G in cluster1 is higher than in cluster2 (Fig. S11A and C). Nine variants per sample were observed in cluster1, and eight variants per sample were observed in cluster2. Compared to cluster2, the mutation rate of primary genes, such as TP53, PIK3CA, RB1, ATM, and CREBBP, is higher in cluster1 (Fig. S11B and D). These results revealed that the patients in cluster1 were subjected higher mutation burden than that in cluster2. We further explored RCD clusters' clinical relevance, mechanism, and immuno-correlation. We observed that the T stage, M stage, histological grade, and pathological stage showed a significantly different distribution in cluster1/cluster2 (Fig. 8D-G). Then the differential expression genes (DEGs) in cluster1 versus cluster2 were calculated and presented to the Clusterprofile R package to estimate the potential mechanisms (Table S18a). Multiple terms of gene ontology enriched in cluster1, such as antigen binding, extracellular matrix, immune receptor activity, inflammatory response, and adaptive immune response (Fig. 8H-J, Table S18b). Multiple immune-related KEGG pathways, including phagosome, systemic lupus erythematosus, complement, and coagulation cascades, were enriched in cluster1 (Fig. 8K, Table S18c). Through the GSEA software analysis, we also found that some hallmark pathways, including inflammatory response, apoptosis, complement, and interferon-gamma response, were highly enriched in cluster1 (Fig. 8L, Table S18d). GSVA also revealed that the activity of these Hallmark pathways was higher in cluster1 than in cluster2 (Fig. S12A). We further analyzed several BLCA-hallmark pathways - urothelial differentiation pathways were enriched in cluster2, whereas the myofibroblasts pathways were enriched in cluster1 (Fig. 8L, Table S18f)48
Immuno-/chemo-/target-therapy Prediction Of Rcd Clusters
Due to multiple immune-related pathways being differently enriched in different RCD clusters, we further the potential immunotherapy response of RCD clusters. We found that the expression of three immune checkpoints (PD1, PDL1, and CTLA4) was higher in cluster1 than in cluster2 (Fig. 9A). Immune and stromal scores were higher in cluster1, whereas tumor purity was lower in cluster1 (Fig. 9B). In addition, the activity of multiple immune cells was higher in cluster1, such as Treg, aDC, and macrophage cells (Fig. 9C, Table S19a). We used the xCell R package further calculate the activity score of the tumor micro-environment, the activity of the immune cell is consistent with Fig. 9C (Fig. S12B, Table S19b). For the stromal cell, the activity of adipocytes, smooth muscle, pericytes, chondrocytes, and fibroblasts was higher in cluster1 than in cluster2 (Fig. 9D, Table S19c). Through the Tumor Immune Dysfunction and Exclusion (TIDE) analysis, we found the TIDE score in cluster1 was lower than in cluster2, and the number of the responder was higher than in cluster2 (P-value < 0.05, Chi-squared test) (Fig. 9E and F, Table S19d). These results all revealed that the patients in cluster1 may respond to the immunotherapy. Furthermore, the estimated IC50 of vinblastine, cisplatin, and paclitaxel in cluster1 was lower than in cluster2, revealing that these patients in cluster1 were sensitive to these drugs (Fig. 9G). We further provided the DEGs between cluster1/cluster2 on the design website, connectivity score closer to 1 indicates that the drug has the most excellent efficacy, and closer to -1 indicates that the drug has the minimal efficacy56. The results showed that doxorubicin (targeting TOP2A) or cytarabine (targeting POLA1, POLB, POLD1, and POLE) were sensitive to the patients in cluster1 (Fig. 9H).
Identification Of Rcd Cluster Prognosis-related Genes
The prognosis difference between RCD clusters made us explore the reason for this difference. So we want to first identify the RCD cluster1/cluster2 prognosis-related genes by conducting the WGCNA analysis. We achieve a scale-free co-expression network by using a power of β = 5 in RCD cluster1 and cluster2 cohorts, respectively (Fig. S13A). The k and p(k) reached or exceeded 0.85, indicating that an optimal scale-free topology network was constructed (Fig. S13B). Then we merge the highly familiar modules by choosing a cut line of 0.25 and minimum module size of 30 (Fig. 11A). Finally, the interest modules of RCD cluster1 (Yellow module: r = 0.49, p = 3e-26) and cluster2 (Brown module: r = 0.62, p = 1e-43; Turquoise module: r = 0.7, 5e-60) were identified (Fig. 11B, Table S21a). Among these, 73 genes were the risk genes, and 47 were the protective genes for MIBC (Fig. 11C, Table S21b and c). Interestingly, among these 73 risk genes, 34 genes come from the yellow module, which is positively correlated with RCD cluster1 (P-value = 1.25e − 05, hypergeometric test)(Fig. 11D). Among the 47 protective genes, 41 genes come from the RCD-cluster2-related module (P-value = 1.26e -05, hypergeometric test) (Fig. 11E). So it is considered that these genes mainly promote the discriminate prognosis between RCD cluster1 and cluster2.
The Relationships Between Rcd Cluster And Consensus Clusters In Bladder Cancer
The previous study has developed six methodologies for the molecular classification of MIBC samples using the mRNA expression data6,57−61. Aure´lie Kamoun, et al. combined these methods and constructed a consensus cluster recently48. To understand the correlation between these molecular classifications and the RCD cluster, we construct these classifications to conduct a comparison (Fig. 11, Table S21). We discovered that RCD cluster1 is highly correlated with basal/squamous clusters, whereas RCD cluster2 is highly correlated with luminal clusters (Fig. 11). The patient with basal subtype usual showed a poor prognosis48, which further explains why the prognosis of the patients in RCD cluster1 is poor. To explore the distribution of these RCD cluster prognosis-related genes, we used a recent single-cell RNA-sequence data published in a recent paper62. We determined 16 clusters based on the data (Fig. 12A). The cell types were annotated by SingleR R-package (Table S22a and b). Among these, most clusters were epithelial cells validated by specific marker-EPCAM, KRT8, and KRT18 (Fig. 12B, Fig. S14A). In addition, the basal or luminal subtype cells were determined by GATA3, KRT5, and CD44 (Fig. S14B)63. We found multiple yellow module risk genes (such as EMP1, ZEB1, F2R, VIM, CAV1, and TIMP1) distributed in fibroblasts (Cluster 13), tissue stem cells (Cluster 13), and endothelial cells (Cluster 16) (Fig. S15). EMP1, ANXA1, SLC7A11, DUSP1, MAP3K5, JUN, TP63, SAR1A, and FTL showed a high expression in basal type cells (Fig. S15). For the RCD cluster2 related module genes, multiple genes showed high expression in luminal type cells, such as P4HB, QPRT, PRDX6, SCD, CISD1, SLC3A2, BAG3, SLC20A1, ERBB3, and so on (Fig. S16 and 17). These results prove the correlation between the RCD cluster and the consensus cluster. In CCLE bladder cancer cells, we also observed a similar phenomenon (Fig. S18). In luminal-type cells, the expression of some RCDGs (such as AIFM3, SRC, ERBB3, ALOX5, ACSF2, PRKCD, BLCAP, and TP63) was higher. In basal-type cells, expression of some RCDGs (such as MLKL, SPHK1, F2R, AKT3, TXNRD1, VIM, and ZEB1) was higher (Fig. 12C, Table S22c).