Molecular subtypes, tumor microenvironment infiltration characterization and prognosis model based on cuproptosis in bladder cancer

Data collection

Figure S1 demonstrated the process of this study.

In Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/), datasets were included based on the following criteria: (1) Histologically diagnosed as BCa; (2) Availability of RNA expression data; (3) Availability of prognostic information; (4) More than 50 tumor samples. In this way, the normalized gene expression data of five BCa datasets (GSE48277, GSE32894, GSE31684, GSE48075 and GSE13507) with relevant prognostic data were obtained from the GEO. The RNA-sequence data (fragments per kilobase million, FPKM) with clinical information of patients were obtained from a BCa dataset in The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov). The gene expression data and corresponding clinicopathological and follow-up information of metastatic urothelial cancer patients who accepted atezolizumab (PD-L1 inhibitor) were obtained from the IMvigor210 trial (http://research-pub.gene.com/IMvigor210CoreBiologies). BCa patients in five GEO datasets and IMvigor210 trial were combined through the “Combat” method from the “sva” R package (https://doi.org/doi:10.18129/B9.bioc.sva) (Leek et al., 2012) to minimize the batch effects. “Combat” method could adjust data for batch effects through an empirical Bayes framework and return an expression matrix that has been corrected for batch effects. As a result, there were 1,306 BCa patients in the combined dataset Table S1 showed the main information on the combined dataset.

Consensus clustering analysis of cuproptosis-related genes

Ten cuproptosis-related genes, including ferredoxin 1 (FDX1), lipoic acid synthetase (LIAS), lipoyltransferase 1 (LIPT1), dihydrolipoamide dehydrogenase (DLD1), dihydrolipoamide S-acetyltransferase (DLAT), pyruvate dehydrogenase E1 subunit alpha 1 (PDHA1), pyruvate dehydrogenase E1 subunit beta (PDHB), metal regulatory transcription factor 1 (MTF1), glutaminase (GLS) and cyclin dependent kinase inhibitor 2A (CDKN2A), were downloaded from a recent study (Tsvetkov et al., 2022). The network of genes biologically associated with cuproptosis-related genes was created in GeneMANIA (http://genemania.org/). Based on cuproptosis-related genes, BCa patients in the combined datasets were classified into distinct molecular subtypes using the R package “ConsensusClusterPlus”. The number of subtypes (K value) was determined according to three criteria: First, the sample sizes of each subtype were relatively equal. Second, the inter-subtype correlation decreased and the intra-subtype correlation increased after clustering. Third, the curve of the cumulative distribution function raised smoothly and gradually.

Associations of molecular subtypes with clinicopathological features

The correlation of molecular subtypes with clinical features and prognoses was investigated. In addition, Kaplan–Meier and log-rank tests were performed using R packages “survminer” and “survival”.

Evaluation of tumor-infiltrating immune cells (TIICs), Tumor immune dysfunction and exclusion (TIDE) and TME scores

The relative proportion of TIICs in each BCa patient was evaluated using CIBERSORT (https://cibersort.stanford.edu/) (Newman et al., 2015). The TIDE values of BCa patients were evaluated on the TIDE website (http://tide.dfci.harvard.edu/) (Fu et al., 2020). The ‘ESTIMATE’ R package was applicated to evaluate the TME scores (https://bioinformatics.mdanderson.org/estimate/index.html) (Yoshihara et al., 2013).

Identification of differentially expressed genes (DEGs)

The DEGs among subtypes were obtained through the R package “limma” with a false discovery rate (FDR) <0.05 and log2—fold-change—>1.5. The potential functions of cuproptosis pattern-related DEGs were explored through the Gene Ontology (GO)/Kyoto Encyclopedia of Genes Genomes (KEGG) analyses in DAVID (https://david.ncifcrf.gov/). Gene Set Enrichment Analysis (GSEA) was performed using GSEA software (http://www.broad.mit.edu/GSEA/) to uncover biological processes based on sets of differentially expressed genes instead of individual genes. Gene Set Variation Analysis (GSVA) was performed using the R package “GSEABase” and “GSVA” to converts gene expression to pathways using gene sets. In GSEA and GSVA analysis, the KEGG (http://www.gsea-msigdb.org/gsea/msigdb/human/genesets.jsp?collection=CP:KEGG) and GO (http://www.gsea-msigdb.org/gsea/msigdb/human/genesets.jsp?collection=GO) gene sets were used. In addition, the cuproptosis score was calculated to quantitatively evaluate the enrichment scores of cuproptosis based on the gene expression of cuproptosis-related genes using the GSVA method.

Cuproptosis gene cluster and cuproptosis-related prognosis signature

We further tried to construct the cuproptosis gene cluster and a cuproptosis-related prognosis signature. First, the univariate Cox regression and Kaplan–Meier analysis were performed on DEGs to obtain the prognostic DEGs. Second, Patients were divided into distinct cuproptosis gene clusters according to the gene expression of the prognostic DEGs using the R package “ConsensusClusterPlus”. Finally prognostic DEGs were introduced into the least absolute shrinkage and selection operator (LASSO) model was used to develop the cuproptosis-related prognosis signature via the R package “glmnet”, and the penalty regularization parameter lambda (λ) was chosen using 10-fold cross validation to obtain genes with non-zero coefficients and minimize the mean square error. Meanwhile, the minimal λ was identified to obtain the genes. The risk scores of the cuproptosis-related prognosis signature were calculated as follows: cuproptosis-related prognosis signature = Coef₁ × expression of gene 1 + Coef₂ × expression of gene 2 + …+ Coef_m × expression of gene m. Coef denotes the corresponding coefficient of the gene.

Chemotherapeutic drugs susceptibility analysis and mutation

The semi-inhibitory concentration (IC50) of BCa patients was evaluated based on the R package “pRRophetic”. The somatic mutations of BCa patients in TCGA were demonstrated by the R package “maftools”.

Development of the nomogram

The cuproptosis-related prognosis signature and clinical factors were subjected to univariate cox analysis. Factors significant in univariate cox analysis were then subjected to multivariate cox analysis to develop a nomogram. Decision curve analysis (DCA), Time-dependent receiver operating characteristic curve (ROC) curves and Calibration plots were applied to explore the performance of the nomogram and clinical factors.

Cell culture and Real-time quantitative PCR (RT-qPCR)

All cell experiments were approved by the Ethics Committee of Shanghai Tenth People’s Hospital (approval number 2021KN108).

BCa cell lines (5637, T24, RT4, J82, and UMUC3) and human normal bladder epithelial cell line SV-HUC-1 were all purchased from the Chinese Academy of Sciences. RT4 was cultured in McCoy’s 5A medium (Thermo Fisher Scientific, Waltham, MA, USA); SV-HUC-1 was cultured in F12k medium (Sigma-Aldrich, St. Louis, MO, USA); J82 was cultured in DMEM (Thermo Fisher Scientific); T24, 5637 and UMUC3 were cultured in RPMI-1640 medium (Thermo Fisher Scientific). All cells were cultured at 37 °C and 5% CO2. All cell culture medias were supplemented with 10% fetal bovine serum (Gibco, Waltham, MA, USA). The protocols of total RNA extraction and RT-qPCR have been described previously (Wang et al., 2021). Data were analyzed with Bio-Rad software (Bio-Rad, Hercules, California, USA). Data were normalized by β-actin expression level in each sample. 2^−ΔΔCt method was used to calculate relative expression changes (Livak & Schmittgen, 2001). Primer sequences are presented in Table S2.

Statistics

Principal component analysis (PCA) was conducted using the R package “pcaMethods”. The correlations among factors were investigated using the R package “corrplot”, “circlize”, “igraph” and “reshape2”. The cutoff points of the cuproptosis score and risk score were determined by the R package “survminer”. SPSS 23.0 (IBM, NY, USA) and R v.4.0.3 software (R Core Team, 2020) and were used for statistical analysis. Variables between subtypes or groups were compared using Mann–Whitney U test, Chi-square test, or Student’s t-test, as appropriate. A P < 0.05 in two-tailed analyses was statistical significance, with a Bonferroni correction applied for pairwise comparisons.

Identification of cuproptosis subtypes in BCa

Figure 1A showed the chromosome locations of the cuproptosis-related genes. GeneMANIA network built around cuproptosis-related genes demonstrated the top 20 genes interacted with cuproptosis-related genes (Fig. 1B). High interactions within cuproptosis-related genes were also showed in this network. The correlation network among cuproptosis-related genes were shown in Fig. 1C. Then, we investigated somatic mutations in cuproptosis-related genes in TCGA and discovered relatively low somatic mutations in these genes (Fig. 1D). Among them, CDKN2A had the highest somatic mutations while others had five or less than five somatic mutations.

Through the “Combat” algorithm, the batch effects were eliminated in the combined dataset (including GSE13507, GSE31684, GSE32894, GSE48075, GSE48277 and IMvigor210 trial) (Figs. S2A, S2B). The gene expression data of cuproptosis-related genes were used to categorize the BCa patients in the combined dataset. As a result, K = 3 was an optimal value to classify the BCa patients in the combined dataset into Cluster 1/2/3 (n = 515/427/364) (Figs. S2C–S2E). PCA plot showed novel differences among the three clusters (Fig. 1E). Survival analysis revealed a poorer overall survival (OS) in Cluster 2/3 patients than those in Cluster 1 patients (P = 0.027; Fig. 1F, Fig. S2F). Ten cuproptosis-related genes had significantly different expression patterns among these three clusters (Fig. 1G, Fig. S2G). In addition, Cluster 2 and Cluster 3 were preferentially related to the high T stage and higher risk of death (Fig. 1H).

Characteristics of cuproptosis score and TIICs infiltration in the cuproptosis subtypes

Survival analysis showed that cuproptosis score was a protective biomarker in BCa patients (Fig. 2A). We then explored the correlation of the cuproptosis score with ten cuproptosis-related genes and cuproptosis subtypes. The results revealed that ten cuproptosis-related genes all had a positive correlation with the cuproptosis score (Fig. 2B, Fig. S3). Cluster 1 patients had the highest cuproptosis scores while Cluster 2 patients had the lowest cuproptosis scores (Fig. 2C).

We further investigated the differences in the proportion of TIICs among cuproptosis subtypes. The relative proportion of the 22 TIICs in each BCa patient was calculated and showed in Figs. 2D, 2E. Eighteen of 22 TIICs had significantly differently composition among cuproptosis subtypes (Fig. 2E), which revealed the diverse tumor immune cell microenvironments among cuproptosis subtypes.

Identification of DEGs and cuproptosis gene cluster

As patients in Cluster 2/3 had similar survival outcomes, we combined Cluster 2 and Cluster 3 (named Cluster 2/3) and identified cuproptosis subtype-related DEGs between Cluster 2/3 and Cluster 1. In this way, 291 cuproptosis subtype-related DEGs were obtained (Fig. 3A), and DEGs that were up-regulated in Cluster 2/3 were enriched in many cancer-promoting biological processes (Figs. 3B, 3C), including Cell division, Cell cycle, DNA biosynthetic process, DNA replication, DNA replication initiation. The cuproptosis subtype-related DEGs were then introduced to univariate Cox analysis and Kaplan–Meier analysis in TCGA and combined dataset to obtain the prognostic DEGs. In this way, 71 prognostic DEGs were the overlapping genes related to OS in TCGA and combined data (Fig. 3D, Table S3). Based on the above genes, BCa patients in the combined dataset were divided into two genomic subtypes using the consensus clustering algorithm, namely GeneCluster 1 and GeneCluster 2 (Figs. S4A–S4F). PCA plot showed significant differences in the cuproptosis transcription profiles among two cuproptosis gene subtypes (Fig. 3E). Patients in GeneCluster 2 had significantly poorer OS than those in GeneCluster 1 (P < 0.001, Fig. 3F), and BCa patients in GeneCluster 2 had significantly higher cuproptosis scores compared with BCa patients in GeneCluster 1 (Fig. 3G). In addition, there were significant differences in clinicopathological features (Fig. S4G), TME characteristics (Figs. S5A, S5B) and cuproptosis-related gene expression (Fig. 3H) between GeneCluster 1 and GeneCluster 2.

Development of the cuproptosis-related prognosis signature

Based on the minimal λ (−2.840), 16 hub genes with their corresponding coefficients were obtained among prognostic DEGs to develop the cuproptosis-related prognosis signature (Fig. 4A, Fig. S6A), in which the partial likelihood deviance was at its minimum (Fig. 4B). The correlations of 16 hub genes in the cuproptosis-related prognosis signature were demonstrated in Fig. S6B, suggesting that these genes were not highly correlated with each other. The risk scores were calculated by summing the genes weighted by corresponding coefficients, which was mentioned above. The optimal cutoff point for risk scores was 0.812 determined by the R package “survminer” in TCGA. BCa patients were classified into high- and low-risk groups according to the optimal cutoff point, and high-risk patients had significantly lower survival rates and poorer OS than low-risk patients (Figs. 4C, 4D). Based on the prognostic DEGs, PCA was performed and revealed that BCa patients in low or high risk groups were distributed in discrete directions indicating novel differences in the prognostic DEGs between the high- and low-risk patients (Fig. 4E). To validate the prognostic value of the cuproptosis-related prognosis signature, BCa patients in the combined dataset were also divided into low- and high-risk groups. Kaplan–Meier curve showed that patients in high risk group also had poor OS than those in the low risk group in the combined dataset (Fig. 4F). The alluvial diagram illustrated the distribution of BCa patients in three cuproptosis subtypes, two risk groups, two gene subtypes, two cuproptosis score groups and survival outcomes (Fig. 4G). Patients in Cluster 2/3 and GeneCluster 1 had higher risk scores compared with those in Cluster 1 and GeneCluster 2, respectively (Figs. 4H, 4I). In TCGA, high-risk patients were preferentially correlated with higher age, higher TNM stage and low survival rate (Fig. 4J, Figs. S6C–S6G).

Characteristics of the signature

We compared somatic mutations, TMB, cuproptosis score and IC50 values of chemotherapeutic drugs between low- and high-risk patients. In TCGA, the top ten mutated genes were demonstrated in Figs. 5A and 5B, respectively. High-risk patients had lower TMB than those in low risk patients (P = 0.03, Fig. 5C). There was a negative correlation between TMB and risk scores (Fig. 5D, R =−0.17, P < 0.001). Similarly, in the combined dataset, low risk patients had higher cuproptosis score compared with high-risk patients (P < 0.001, Fig. 5E), and risk score was negatively related to the cuproptosis score (Fig. 5F). In addition, among seven common chemotherapeutic drugs used in BCa, the IC50 of five drugs were significantly different between two risk groups (Fig. 5G). GSEA revealed that several cancer-promoting pathways were highly enriched in high risk patient, including ‘Focal adhesion’, ‘MAPK signaling pathway’ and ‘Pathway in cancer’, while the ‘Positive regulation of necrotic cell death’ was enriched in low risk BCa patients (Fig. 5H). GSVA showed similar results with the ‘MAPK signaling pathway’, ‘JAK STAT signaling pathway’, ‘Focal adhesion’ and ‘Pathway in cancer’ and ‘WNT signaling pathway’ enriched in high risk BCa patients (Fig. 5I).

Characteristics of the TME in distinct risk groups

Fourteen of 22 TIICs had different compositions between two risk groups (Fig. 6A, Table S4), which revealed the diverse tumor immune cell microenvironments between risk groups. Specifically, the proportions of macrophages were higher in high-risk patients compared with low risk patients, while T cells follicular helper, T cells CD8, B cells memory and monocytes contained significantly lower infiltration levels in high-risk patients compared to those in low risk patients. In addition, the risk score had a significant correlation with 12 TIICs (Fig. 6B, Fig. S7). We then investigated the relationship between the expression of some immune checkpoints and the risk score. The results revealed that some immune checkpoints had different expression levels between two risk groups, including CD274, PDCD1 and CTLA4 (Fig. 6C), and risk score had a correlation with some immune checkpoints (Fig. 6D). High-risk patients had higher immune scores, stromal scores and TIDE scores compared with low-risk patients (Figs. 6E, 6F). Based on the TIDE scores, high-risk patients had a lower proportion of responders compared to those in low-risk patients if received immunotherapy (Figs. 6G, 6H). These results indicated that this model could be a useful predictor of the efficacy of immunotherapy. GSEA and GSVA revealed that high-risk patients had multiple immune-related pathways enriched, including ‘Regulation of T cell activation’, ‘Negative regulation of immune system process’, ‘Regulation of T cell differentiation’, ‘T cell receptor signaling pathway’ and ‘TGF beta signaling pathway’ (Figs. 5I, 6I).

Construction of the nomogram

To overcome the inconvenient clinical utility of the signature (risk score), a cuproptosis-related prognosis signature-based nomogram was constructed to predict the survival outcomes of BCa patients. The results of univariate and multivariate Cox regression analysis showed that risk score and T stage remained the independent prognostic factors (Fig. 7A), so the signature and T stage were selected to develop the nomogram (Fig. 7B). The histogram also showed that patients with T stage ≥T3 and high-risk patients had higher death rates than others (P < 0.001, Fig. 7C). The nomogram achieved the optimal performance in the prediction of patients’ survival outcomes compared with other factors (Figs. 7D–7H). Furthermore, the nomogram achieved higher clinical net benefit than the T stage and risk group (Figs. 7I–7K).

The experimental validation of the cuproptosis-related genes

The mRNA expression pattern of the ten cuproptosis-related genes in five human BCa cell lines (T24, 5637, RT4, J82, and UMUC3) and human normal bladder epithelial cell line SV-HUC-1 was determined using RT-qPCR. The results showed significant differences in the expression levels of nine cuproptosis-related genes (except LIAS) between BCa and normal bladder epithelial cells (Fig. 8).