ADME gene-driven prognostic model for bladder cancer: a breakthrough in predicting survival and personalized treatment

Abstract

Background

Genes that participate in the absorption, distribution, metabolism, excretion (ADME) processes occupy a central role in pharmacokinetics. Meanwhile, variability in clinical outcomes and responses to treatment is notable in bladder cancer (BLCA).

Methods

Our study utilized expansive datasets from TCGA and the GEO to explore prognostic factors in bladder cancer. Utilizing both univariate Cox regression and the lasso regression techniques, we identified ADME genes critical for patient outcomes. Utilizing genes identified in our study, a model for assessing risk was constructed. The evaluation of this model's predictive precision was conducted using Kaplan–Meier survival curves and assessments based on ROC curves. Furthermore, we devised a predictive nomogram, offering a straightforward visualization of crucial prognostic indicators. To explore the potential factors mediating the differences in outcomes between high and low risk groups, we performed comprehensive analyses including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG)-based enrichment analyses, immune infiltration variations, somatic mutation landscapes, and pharmacological sensitivity response assessment etc. Immediately following this, we selected core genes based on the PPI network and explored the prognostic potential of the core genes as well as immune modulation, and pathway activation. And the differential expression was verified by immunohistochemistry and qRT-PCR. Finally we explored the potential of the core genes as pan-cancer biomarkers.

Results

Our efforts culminated in the establishment of a validated 17-gene ADME-centered risk prediction model, displaying remarkable predictive accuracy for BLCA prognosis. Through separate cox regression analyses, the importance of the model’s risk score in forecasting BLCA outcomes was substantiated. Furthermore, a novel nomogram incorporating clinical variables alongside the risk score was introduced. Comprehensive studies established a strong correlation between the risk score and several key indicators: patterns of immune cell infiltration, reactions to immunotherapy, landscape of somatic mutation and profiles of drug sensitivity. We screened the core prognostic gene CYP2C8, explored its role in tumor bioregulation and validated its upregulated expression in bladder cancer. Furthermore, we found that it can serve as a reliable biomarker for pan-cancer.

Conclusion

The risk assessment model formulated in our research stands as a formidable instrument for forecasting BLCA prognosis, while also providing insights into the disease's progression mechanisms and guiding clinical decision-making strategies.

Introduction

Annually, bladder cancer (BLCA) accounts for about 3% of new global cancer diagnoses, ranking it as the 11th most common cancer. It is estimated that each year sees around 573,000 new diagnoses [1]. BLCA arises from a complex interplay of factors such as genetics, environmental and occupational influences, tobacco use, and conditions linked to obesity [2]. In terms of genetics, the notorious Lynch syndrome increases the risk of developing BLCA [3]. Smoking, especially black smoke, leads to BLCA mainly because of the high amount of carcinogens in it and is metabolized by the kidneys and enriched in the bladder [4]. As for occupational exposure, it has been reported that the rate of BLCA is significantly higher in workers who are chronically exposed to aromatic amines, PAH, such as petroleum workers, metal factory workers [5]. Up to 2024, the 5-year survival rate of BLCA is encouraging, nearing 78%, suggesting a favorable prognosis, but is the fourth leading cause of cancer deaths in men over 80 years of age [6]. Nevertheless, metastatic BLCA presents a sharply reduced overall survival (OS) rate, which dips under 6% [7]. Significantly, around half of those initially diagnosed with muscle-invasive BLCA manifest recurrences and progress to metastasis [8]. In addition, the extensive exploration of cancer biology since the past centuries has facilitated the implementation of precision treatments, such as integrating surgery and chemotherapy, selecting appropriate molecular drugs, and breaking through the barriers of the immune microenvironment [9]. In summary, the malignant evolution of bladder cancer and the progress in precision cancer therapy emphasize the urgent need to construct reliable biomarker models.

The drug ADME (absorption, distribution, metabolism, excretion) gene superfamily is instrumental in the intricate mechanisms of drug absorption, distribution, metabolism, and elimination [10,11,12]. Per the PharmaADME consortium (http://www.pharmaadme.org), the superfamily encompasses 266 genes, including a fundamental subset of 32 genes classified by their specific functions in pharmacokinetics including drug transporters, regulatory factors, and enzymes active in phase I and II drug metabolism [13,14,15]. In addition to affecting individual drug response [16], polymorphisms of ADME genes are closely associated with many aspects of cancer. For example, single nucleotide polymorphisms in the metabolic gene AKR1C3 have been associated with bladder cancer [17] and leukemia [18], polymorphisms in genes such as CYP17A1 affect drug efficacy in prostate cancer [19], and 3'UTR alterations in a group of ADME genes promote the risk of toxicity in breast cancer treatment [20]. Moreover, significant intergroup variations exist in the transcriptional, translational, and epigenetic modulation of ADME gene expression [21,22,23].

Recent studies have provided insights into the prognostic potential of ADME heterogeneity in renal cancer from the perspective of spatial transcriptome combined with single-cell data as well as bulk RNA-seq, deepening the understanding of the regulatory role of ADME genes in the tumor microenvironment [24, 25]. Hu and associates have pinpointed critical ADME genes that predict overall survival in a range of cancers [26]. Research by Han et al. revealed connections between ADME gene variations and survival outcomes in non-small cell lung carcinoma [27]. Tang and his team have associated these genes with prognostic results and the extent of immune cell infiltration in head and neck squamous cell carcinoma [28].

Although progress has been made in understanding the function of ADME genes in cancer, the potential clinical significance of genetic anomalies related to ADME in BLCA has not been fully explored. In this research endeavor, we sought to unravel the molecular alterations and clinical significance of ADME-associated genes in BLCA. We introduced a novel risk prediction model to evaluate treatment outcomes and prognosis in BLCA patients. Independent cohorts were used to rigorously validate this model, confirming its strong predictive capabilities and its potential to influence clinical decisions. The deployment of this model could imminently facilitate the development of customized treatment strategies for patients diagnosed with BLCA.

Materials and methods

Data collection

Clinical information and transcriptional profiles from the TCGA-BLCA and TCGA-Pan-cancer dataset were retrieved through the UCSC Xena platform (https://xena.ucsc.edu) [29]. During our analysis, any subjects missing complete survival data or RNA sequencing details were systematically excluded. We obtained the GSE13507 gene expression dataset [30] as independent validation set, which includes 165 primary bladder cancer specimens, from the Gene Expression Omnibus (GEO) repository. The GSE13507 cohort was widely used as a validation set for BLCA prognostic models [31,32,33]. Data from the IMvigor 210 cohort, detailing clinical and characteristic information, is accessible through the 'IMvigor210CoreBiologies' R package. 298 ADME genes were obtained from the supplementary table in the previous literature of Tang et al. [28].

Development and validation of a prognostic risk prediction model using ADME genes

In our study of the TCGA-BLCA cohort, based on the threshold selection strategy in previous literature [24, 34], by setting a log2FoldChange threshold greater than 1 and an adjusted P-value under 0.05, researchers identified genes exhibiting differential expression between tumor tissues and adjacent normal tissues through “edgeR” R package. To evaluate their prognostic significance, a univariate Cox regression analysis was performed on ADME-related differentially expressed genes (DEGs) [34]. Lasso regression is known for its efficiency in screening important features and can effectively deal with the problem of multicollinearity of parameters, avoiding overfitting to a great extent [35]. In our study, a robust risk assessment model was developed by identifying the optimal regularization parameter, lambda (λ), through Lasso cox regression analysis paired with tenfold cross-validation. This approach ensured precise determination of λ, enhancing the model’s efficacy. Risk scores for patients were derived by integrating standardized gene expression values (Expi) with their respective regression coefficients (Coei), following a strictly defined algorithmic framework:

$$Risk score=\sum_{\text{i}=1}^{\text{n}}(Expi \times Coei)$$

Risk categorization of samples into high-risk and low-risk groups was performed using the median risk score as the dividing criterion. Drawing on the analytic strategies of previous studies, the disparities in survival between the categories were assessed through Kaplan–Meier survival curves and analyzed using log-rank tests [36,37,38]. To confirm the predictive accuracy of our risk model, time-dependent ROC curves were constructed [39, 40]. In addition, Decision Curves Analysis (DCA) was performed to aid in the understanding of the prognostic potential of the model [38]. The robustness of our model's predictive capability was further validated using BLCA samples sourced from the GEO database, emphasizing its consistent performance across diverse datasets, thereby enhancing its credibility and dependability.

Establishment and validation of the nomogram

Cox proportional hazards regression analysis was utilized to ascertain the independent prognostic capability of the risk score. A nomogram was formulated, integrating key clinical parameters—age, gender, pathological stage, and the risk score—to predict outcomes, aimed at refining the precision of prognostic predictions [40, 41]. The predictive accuracy of this nomogram was rigorously scrutinized using calibration curves, DCA curves and ROC analysis, underscoring its reliability and clinical applicability in predicting patient outcomes.

Enrichment analysis in high- and low-risk groups

By setting a log2FoldChange threshold above 1 and keeping the adjusted P-value below 0.05, we pinpointed genes that were differentially expressed between high- and low-risk groups. These DEGs underwent extensive functional characterization using the GO and KEGG databases to unravel their biological significance. Furthermore, we performed Gene Set Enrichment Analysis (GSEA) tailored to GO and KEGG datasets, highlighting unique metabolic pathways enriched in each risk group, providing deeper insights into the molecular underpinnings of their distinct risk profiles. The results of enrichment analysis are presented according to adjusted P-value.

Immune infiltration analysis

Using the 'ESTIMATE' R package, we derived stromal, immune, and overall ESTIMATE scores to delineate the cellular composition of the tumor microenvironment (TME) in each specimen. To characterize the TME in BLCA patients, using the CIBERSORT algorithm [42] in conjunction with ssGSEA methods, the extent of immune cell infiltration was quantified. Following quantification, Spearman correlation analysis was utilized to investigate the association between immune cell infiltration and the risk score. In addition, we analyzed the differences in the expression of human leukocyte antigen (HLA) family members and immune function assessed by ssGSEA between high and low risk subgroups.Finally we uncovered correlations between immune checkpoint expression and risk scores and explored the potential of risk scores as a prognostic indicator for immunotherapy in the cohort Imvigor210 treated with anti-PD-L1.

Malignancy analysis, somatic mutation analysis and evaluation of chemotherapeutic drug sensitivity

By ssGSEA, we assessed the Angiogenesis score as well as the epithelial-mesenchymal transition (EMT) score for each sample. By applying the “maftools” R package, we carefully analyzed the data on somatic mutations and compared tumor mutation burden (TMB) in different risk subgroups. KM curves were used to demonstrate survival differences between subgroups. We utilized the 'oncoPredict' R package to evaluate differential drug responsiveness across high- and low-risk cohorts, revealing potential therapeutic approaches informed by our findings.

Construction of ADME-related clusters

“ConsensusClusterPlus” R package was utilized to identify diverse ADME patterns and cluster patients in TCGA cohort into multiple categories for further study. Principal component analysis (PCA) and tSNE were applied to distinguish clustered differences in expression profiles. KM curves were drawn to demonstrate the prognostic characteristics of different clusters. In addition, we explored differences in pathway activation across clusters by GSEA and GSVA.

Identification of core gene

Utilizing the STRING database (http://string-db.org), we constructed a PPI network from genes we identified as prognostic markers.The protein that communicates most extensively with other proteins was identified as the core of the network. Core gene-associated ppi networks were constructed with GeneMANIA (https://genemania.org/). We assessed the prognostic potential of the core genes by KM curves and revealed possible means of regulation by immune infiltration assessment as well as GSVA analysis. In addition, we screened for differential genes between high and low expression groups by the same criteria as before, and based on this we performed KEGG enrichment analysis to deepen our understanding of the mechanism.

Cell culture

Four BLCA cell lines (T24, UMUC3, RT4, and BIU87) and normal cells (SVHUC) were obtained from the Cell Bank of the Chinese Academy of Sciences (Shanghai, China). BIU87 cells were treated with 10% PMI-1640 medium, T24 cells with 10%McCoy's 5A medium, RT4 cells with 10% DMEM medium, and UM-UC-3 cells with 10% MEM medium. SV-HUC-1 was cultured using 10% Ham's F-12 K medium in a constant temperature incubator at 37 °C with 5% CO2 volume fraction using culture bottles. All were cultured at 37 °C in a humidified atmosphere containing 5% CO₂.

Quantitative real-time PCR (qRT-PCR) and immunohistochemical analysis

Our previous study described specific methods for total RNA extraction and qRT-PCR. The following primers were used in this study:

• CYP2C8 forward, 5'-CATTACTGACTTCCGTGCTACAT-3' and

• CYP2C8 reverse, 5'-CTCCTGCACAAATTCGTTTTCC-3';

• Actin forward, 5'-GAAGATCAAGATCATTGCTCCTC-3' and

• Actin reverse, 5'-ATCCACATCTGCTGGAAGG-3’.

Actin was used as an internal control. The procedures were performed three times to ensure accuracy and precision. The relative expression levels were calculated using the 2^−ΔΔCT method. Immunohistochemical images for the gene were sourced from The Human Protein Atlas, comparing both cancerous and non-cancerous bladder tissues.

Core gene-based pan-cancer analysis

Pan-cancer differential expression plot of the core genes were downloaded from Timer website. Then by univariate cox analysis, we analyzed the potential of the core genes as predictors of OS, Progression Free Interval (PFI), and Disease Free Interval (DFI) in the TCGA-Pan-cancer dataset, and demonstrated the KM curves in a number of significantly discriminatory cancer types. Subsequently, we explored the correlation of core gene expression with immune checkpoint expression and abundance scores derived from the Cibersort algorithm.In addition, we explored the correlation between core gene expression and scores based on KEGG's immune-inflammation-related gene set after assessing the enrichment score by the ssGSEA method.

Statistical analysis

To assess predictive accuracy, the analysis was conducted using ROC curves, utilizing the AUC as the key quantitative measure. Additionally, we employed the Wilcoxon test for nonparametric analysis. In analyzing the data, the software R, version 4.3.3, was employed. The significance of differences among groups was determined with a threshold P-value set at below 0.05, while the significance test for DEGs was then based on an adjusted P-value less than 0.05. P-value correction was done by the Benjamini-Hochberg (BH) method.

Results

Screening of differentially expressed and prognosis-associated ADME genes in BLCA

The specific steps of the study are detailed in the flowchart (Fig. 1). Utilizing the RNA-seq data from the TCGA dataset, we meticulously analyzed the transcriptional profiles of BLCA and corresponding normal samples, uncovering a total of 6444 genes that displayed significant dysregulation (adjusted P-value < 0.05), including 80 ADME genes. The prognostic significance of specific ADME genes was assessed using univariate Cox regression analysis. This evaluation demonstrated that 19 of these genes significantly correlate with patient outcomes (Fig. 2 A). We applied Lasso regression to enhance the precision of the analysis, effectively reducing the list of prognostic ADME genes from 19 to 17 key genes associated with prognosis: CYP2C8, SLC22A1, ABCA4, ABCC9, CHST6, CYP1B1, CYP27B1, CYP3A43, CYP3A7, CYP4F8, NR1I2, PON1, SLC22A11, SLC22A3, SULT1A3, TAP1, and UGT2B4 (Fig. 2 B-C).

Fig. 1

Flowchart of this study

Fig. 2

Construction of ADME-related BLCA risk model. A Forest plot capturing 19 ADME genes pinpointed through univariate Cox regression as potential prognostic indicators. B and C Lasso regression analysis reduces dimensions, focusing on these 19 genes. D and I Categorizing both the training cohort and the validation cohort into high-risk and low-risk groups based on risk scoring. E and J Scatter plots comparing survival time to survival outcomes within the training and validation cohorts. F and K Kaplan–Meier curves depicting OS for high- and low-risk groups in both cohorts. G and L The ROC curves for prognosis prediction based on risk scores in both cohorts. H and M Heatmaps illustrating the expression levels of prognosis-related ADME genes within high-risk and low-risk groups for both cohorts

Development and validation of a prognosis-related ADME gene risk prediction model

Utilizing multivariate Cox analysis, we constructed predictive models and obtained risk scores for each sample. In the TCGA training dataset for BLCA patients, we divided the subjects into high-risk and low-risk groups based on the median risk score (Fig. 2 D). The scatterplot visualized the increase in the number of dead patients and the decrease in survival time as the risk score increased (Fig. 2 E), while the Kaplan–Meier curves highlighted the significant OS differences between the high and low risk groups (Fig. 2 F). This certainly demonstrated the excellent prognostic stratification ability of the model. The ROC curves were used to evaluate our model's predictive accuracy (Fig. 2 G), revealing AUC scores of 0.697, 0.681, and 0.713 for the 1-, 3-, and 5-year overall survival forecasts, respectively. We visually details the expression patterns of 17 ADME genes associated with prognosis, displayed across two risk clusters in a heatmap format (Fig. 2 H). To fortify the credibility of our risk prediction model, we executed validation using an independent GEO dataset (GSE13507, n = 165). After distinguishing between high and low risk groups (Fig. 2 I), scatter plots (Fig. 2 J) alongside Kaplan–Meier curves (Fig. 2 K) from the test set illustrated significantly shorter overall survival times for patients in the high-risk group when compared to their low-risk counterparts. We then depicts the validation cohort's predictive outcomes, showcasing the AUC scores of 0.677, 0.669 and 0.695 for the 1-, 3-, and 5-year overall survival predictions, respectively (Fig. 2 L). Furthermore, to delineate the expression variations of the 17 ADME genes linked to prognosis across different risk levels, a heatmap was utilized (Fig. 2 M). In addition, we plotted DCA curves and calibration curves in the TCGA as well as the GEO cohort (Fig S1), and the high sensitivity further reflected the potential of the model to guide the process of clinical application.

Identification of clinical features of the ADME-related model

Differences in the distribution of clinical characteristics in high- and low-risk subgroups were exquisitely illustrated in pie charts and heatmap (Fig. 3A-B). Initially, we compared risk scores in cohorts of patients with different clinical characteristics. It was found that although risk scores did not differ significantly between male and female genders, patients of advanced age and those with advanced pathologic staging demonstrated more significantly higher risk scores (Fig. 3 C). In addition, we analyzed the distribution of clinical characteristics in different risk subgroups (Fig. 3 D). Interestingly, there was a significantly higher proportion of patients with advanced pathologic stage as well as advanced age in the high-risk group. Subsequent KM survival analyses demonstrated the power of ADME-related model to differentiate patient survival across diverse cohorts, and we found that highly scored patients survived worse across patients with different clinical characteristics (Fig. 3 E). These results indicate a significant negative correlation between ADME-related score and prognosis in BLCA patients, which further validated the prognostic accuracy of our model.

Fig. 3

Association between ADME-related signature and clinical features. A and B Pie chart and heatmap showing the distribution of people with different clinical characteristics in the high and low risk groups. C Comparison of ADME-related score in different clinical groups. D Proportional differences in clinical characteristics between high- and low-risk groups. E Validating the prognostic accuracy of our model in cohorts with different clinical characteristics

Prognostic significance of the risk score

To assess the risk score's effectiveness as an independent predictor for BLCA, both univariate and multivariate Cox regression analyses were performed. Statistical analysis revealed that, within the training cohort, the risk score significantly influenced prognosis (Fig. 4 A and B). This finding was subsequently corroborated in an external cohort, reinforcing the reliability of the risk score as a prognostic biomarker (Fig. 4 C and D). To refine prognostic precision, we merged the training set with the validation set to construct a nomogram, incorporating the risk score with variables such as gender, age, and stage (Fig. 4 E). The model's predictive precision was demonstrated by a calibration plot that matched the ideal trajectory and DCA curves (Fig. 4 F and G). The evaluation of the nomogram through ROC analysis produced AUC values of 0.795, 0.788, and 0.799 for predictions of 1-, 3-, and 5-year overall survival, respectively (Fig. 4 H), which underscores the nomogram’s robust predictive capability.

Fig. 4

Risk score as an independent prognostic indicator for BLCA patients, alongside the development and confirmation of a Nomogram. A and C Forest plots from univariate Cox regression analysis within the TCGA and GEO cohorts. B and D Forest plots from multivariate Cox regression analysis for both cohorts. E A Nomogram designed to forecast OS rates in BLCA patients. F Calibration curves associated with the Nomogram. G DCA curves associated with Nomogram features. H The ROC curves evaluating the Nomogram's predictive performance

Functional enrichment analysis

In-depth investigations leveraged GO and KEGG databases to elucidate the functional distinctions of genes with varying expression profiles across high- and low-risk cohorts. GO analysis identified significant clustering in biological functions, specifically concerning the organization of the extracellular matrix and associated functions. Key areas include the architecture of collagen-rich extracellular matrices, activation of signaling receptors, and interactions between receptors and ligands (Fig. 5 A-B). Moreover, the KEGG analysis illuminated the enrichment of DEGs within a series of cancer-promoting pathways including PI3K − Akt signaling pathway, Proteoglycans in cancer and Focal adhesion (Fig. 5 C). To delve deeply into the functional pathways, GSEA was also utilized. In this analysis, immune-related pathways, particularly those involving response to type-I interferon, antigen processing and presentation, toll-like receptor signaling cascades showed a significantly higher prevalence in the low-risk group (Fig. 5 D). We observed significant enrichment of pathways concerning immune function, notably leukocyte migration and immunoglobulin complex activities, in the high-risk group (Fig. 5 E). These distinctions outline unique immune signatures for the two groups.

Fig. 5

Enrichment analyses for high- and low-risk categories. A and B The GO enrichment assessments of DEGs. C The KEGG analysis. D and E GSEA between the two groups

Immune infiltration analysis

The above studies established a robust link between the ADME-associated risk score and immunological pathways, we therefore performed careful analyses of the tumor immune microenvironment. The high-risk subset was characterized by increased stromal, immune, and ESTIMATE scores (Fig. 6 A). Utilizing the cibersort algorithm enabled the evaluation of diverse immune cell subsets in the TCGA training cohort, enhancing the understanding of the tumor's immune landscape. We presented a boxplot that displayed the infiltration levels of multiple immune cells (Fig. 6 B). Notably, there were elevated concentrations of CD8 + T lymphocytes, in the low-risk group, the higher presence of activated dendritic cells and ollicular helper T cells were noted as well. Conversely, the high-risk group exhibited increased levels of naive B cells, quiescent CD4 + memory T cells, and both M0 and M2 macrophages, as well as resting mast cells. Utilizing the ssGSEA approach, the prevalence of various immune cell types was quantified and visually represented through a heatmap (Fig. 6 C). Spearman's correlation analyses emphasized significant linkages between the risk score and patterns of immune cell infiltration (Fig. 6 D). Moreover, we present the correlation between model genes and the abundance of immune cell infiltration in Fig S2. Interestingly, we found that the expression of many important HLAs was down-regulated in high-risk subgroups, such as HLA-A, B, and C (Fig. 6 E). Subsequently, we revealed differences in immune function between high- and low-risk subgroups, with MHC I-associated functions being more active in the low-risk subgroup and macrophage-associated functions being significantly up-regulated in the high-risk subgroup (Fig. 6 F). Moreover, we illustrates a strong inverse correlation between the risk score and the expression levels of various immune checkpoint genes, such as CD40, TNFRSF14 and TNFRSF25 (Fig. 6 G). The discovery of distinctive infiltrations of immune cells in subsets with varying risk levels led to a deeper analysis of the IMvigor 210 cohort, focusing on exploring the relationship between ADME scores and the efficacy of immunotherapy treatments. Findings indicated that individuals in the low-risk subset demonstrated not only better prognosis (Fig. 6 H), but also enhanced responses to immunotherapy (Fig. 6 I), suggesting a higher likelihood of benefit from immunotherapy in this subgroup.

Fig. 6

Immune infiltration analyses and immunotherapeutic analyses. A The TME scores contrasting low and high-risk patients. B Boxplots of the infiltration levels of 22 specific immune cells. C Heatmap visually differentiating the quantities of 28 immune cells between high- and low-risk groups. D The relationships between risk scores and the infiltration of immune cells. E Differential expression of HLA-related genes. F Landscape of immune function differences. G The correlation between risk score and immune checkpoints, assessed using the Pearson test, where a significance of P < 0.001 was recorded. H The Kaplan–Meier curve comparing OS between the two risk subgroups in IMvigor210 cohort. I Box plots examining the relationship between risk score and the effectiveness of immunotherapy

Malignancy as well as mutation analysis

Initially, we analyzed the angiogenesis scores obtained based on ssGSEA as well as the EMT scores. The results showed that both were significantly upregulated in the high-risk group and had an extremely significant positive correlation with the risk score (Fig.7 A-D). Subsequently, we assessed differences in TMB across risk groups. Interestingly, we found that high TMB was enriched in the low risk group (Fig. 7 E) and revealed a negative correlation between TMB and risk scores (Fig. 7 F). Further KM analysis demonstrated good prognostic discriminatory validity of the TMB (Fig. 7 G), leading us to analyze the TMB in conjunction with ADME-related score. As we expected, the joint analysis demonstrated strong prognostic discrimination, indicating that patients with high risk and low TMB exhibited the worst prognosis, while those with low risk and high TMB showed the best prognosis (Fig. 7 H). In addition, we present the mutational landscape of the whole sample (Fig. 7 I). Moreover, we present not only the mutations of the 17 prognostic ADME genes (Fig. 7 J), but also the mutation ranking of the genes in the high and low subgroups in exquisite waterfall plots (Fig. 7 K-L). Interestingly, the pathway with the highest mutation rate was statistically found to be consistent in both the high and low risk subgroups, both being the RTK-RAS pathway. In addition, the low-risk group had a higher mutation rate of the WNT pathway, whereas the mutation rate of the Notch pathway in the high-risk group was much higher than that in the low-risk group (Fig. S3).

Fig. 7

Malignancy and mutation analyses. A The difference of angiogenesis score in the two risk subgroups. B Correlation between angiogenesis score and risk score. C The difference of EMT score in the two risk subgroups. D Correlation between EMT score and risk score. E The difference of TMB in the two risk subgroups. F Correlation between TMB and risk score. G Survival curves of the two TMB ssubgroups. H Survival analysis between different subgroups classified by ADME-related risk score and TMB. I Overall mutation landscape of BLCA patients. J Waterfall plot demonstrates the mutation landscape of 17 prognostic genes in tumor samples. K-L Waterfall plots show the somatic mutations in the high- and low-risk subgroups, respectively

Evaluation of medication reactions

Between the cohorts classified as high- and low-risk, we then examined the variations in drug response. As illustrated in Fig. 8, the high-risk subset demonstrated reduced IC50 values for Staurosporine, Obatoclax Mesylate, Foretinib, Dasatinib, and AZ960 compared to the low-risk subset, suggesting an enhanced sensitivity to these medications. In contrast, the low-risk cohort displayed markedly reduced IC50 values for therapeutic agents such as Gefitinib, Afatinib, Oxaliplatin, Palbociclib, Leflunomide, and Cyclophosphamide, suggesting increased sensitivity to these drugs within this group. Comprehending the sensitivity profiles of both immunotherapeutics and chemotherapeutics provides crucial information that can facilitate the customization of treatment strategies tailored to individual patient characteristics.

Fig. 8

Box plots contrasts the responses to chemotherapy across two patient cohorts categorized by risk score

Cluster analysis based on prognostic genes

Based on the expression of the 17 ADME genes, we clustered the TCGA samples into two categories for further exploration of the guidance potential of prognostic genes (Fig. 9 A). Plots of PCA and tSNE demonstrates that clustering has a good discriminatory effect on expression profiles (Fig. 9 B-C). Interestingly, the sankey diagram demonstrated that the C1 subgroup was pointed more toward high risk (Fig. 9 D). We then found that the C1 cluster exhibited higher risk scores (Fig. 9 E–F) and had a significantly lower prognosis than the other cluster (Fig. 9 G), corroborating the prognostic potential of the 17 ADME genes. In view of this, we tried to explore the pathway activation in different clusters. KEGG-based GSEA showed that pathways such as focal adhesion and cytokine-cytokine receptor interaction were significantly enriched in the C1 cluster, whereas the C2 cluster showed more intense metabolic dysregulation (Fig. 9 H-I). In addition, Hallmark-based GSVA heatmaps demonstrated significant upregulation of various pro-oncogenic pathways in the C1 cluster (Fig. 9 J).

Fig. 9

Analyses of ADME-related clusters. A Different ADME-related clusters were identified for k = 2 in TCGA cohort. B PCA analysis of the two clusters. C tSNE analysis of the two clusters (D) Sankey plot of clusters, risk and outcome. E The difference of risk score in two clusters. ****P < 0.0001. F Proportional differences of risk in different clusters. G Overall survival difference between clusters. H-I GSEA analysis in different clusters. J Heatmap of GSVA

Identification and analysis of core gene

Firstly a PPI network was constructed from String database to investigate the functional relationships among 17 ADME genes associated with prognosis. This network illustrated complex interrelationships among the proteins involved (Fig. 10 A). Notably, CYP2C8 emerged as a central node within this network, suggesting its pivotal role. We then constructed a CYP2C8-centered PPI network through GeneMANIA, and interestingly, EPHX2 was the protein with which it most closely interacted (Fig. 10 B). Subsequent survival analyses showed CYP2C8 to be a reliable predictor of OS (Fig. 10 C) and PFI (Fig. 10 D) in the TCGA-BLCA cohort and were validated in the GSE13507 cohort (Fig. 10 E), where we found a better prognosis in the high expression group. Furthermore, we demonstrated different degrees of immune cell infiltration in the high and low CYP2C8 expression groups (Fig. 10 F) and analyzed the correlation between expression and infiltration abundance (Fig. 10 G). Interestingly, activated dendritic cells showed an extremely strong positive correlation with CYP2C8 expression. In addition, the heatmap of GSVA demonstrated significantly upregulated pro-oncogenic pathways in the low-expression group, which coincided with its poorer prognosis (Fig. 10 H). Subsequent KEGG enrichment based on differential genes between the high and low expression groups indicated that cytokine-related pathways as well as metabolic pathways were significantly dysregulated between the two groups (Fig. 10 I).

Fig. 10

Identification and analysis of core gene (A) PPI network of prognostic genes (B) PPI network around CYP2C8 (C) OS segmentation of the TCGA cohort based on CYP2C8 expression (D) DFI segmentation of the TCGA cohort based on CYP2C8 expression (E) OS segmentation of the GSE13507 cohort based on CYP2C8 expression (F) Different immune landscapes in high and low expression groups (G) Correlation of CYP2C8 expression with immune cell abundance (H) Heatmap of GSVA in high and low expression groups (I) KEGG enrichment of DEGs between high and low expression groups

In the previous differential gene expression analysis, CYP2C8 was extracted as a gene whose expression was upregulated in bladder cancer. Utilizing immunohistochemistry images from the HPA database, our analysis focused on comparing BLCA tissue to normal bladder tissue to confirm the result at the protein level (Fig. 11 A). Furthermore, encouragingly, qRT-PCR of cell lines as well as tissue samples verified higher expression of CYP2C8 in bladder cancer (Fig. 11 B-C). These analyses consistently affirmed the enhanced expression of CYP2C8 in bladder cancer samples, further emphasizing its potential importance in this disease.

Fig. 11

Expression validation of CYP2C8 (A) The immunohistochemistry image depicting the expression of CYP2C8 in normal bladder tissue and BLCA tissue (B) Differences in mRNA expression of CYP2C8 in normal and BLCA cell lines (C) Differences in mRNA expression of CYP2C8 in normal and BLCA tissues

Core gene-based pan-cancer analysis

We searched a large number of literatures and found that few studies have investigated CYP2C8 as a tumor biomarker, so we attempted to perform a pan-cancer analysis of it. Initially, we showed the pan-cancer expression landscape of CYP2C8 (Fig. 12 A). Then by univariate cox analysis, CYP2C8 was determined to be significantly associated with OS, DFI and PFI in a range of cancers, such as PAAD, KIRP, GBMLGG, etc. (Fig. 12 B-D). Subsequently, KM curves were utilized to more visually demonstrate the prognostic potential of CYP2C8 in a number of cancer types with significant survival differences between high and low expression groups (Fig. 12 E–G). In addition, we shed some light on the potential regulatory direction of CYP2C8. We accordingly analyzed the correlation of CYP2C8 expression with immune cell infiltration (Fig. 13 A), immune checkpoint expression (Fig. 13 B), immune inflammation-related functions (Fig. 13 C) and stemness score (Fig. 13 D). We have a number of important findings, such as the fact that almost all immune checkpoint genes in KIRC showed a positive correlation with CYP2C8, whereas the opposite was true in KICH. Finally, we comparatively verified the differential expression of CYP2C8 in some cancer types by immunohistochemical profiles derived from the HPA database (Fig. 13 E).

Fig. 12

Prognostic potential of CYP2C8 in pan-cancer. A Differential pan-cancer expression of CYP2C8 (B-D) Prognostic potential of CYP2C8 for OS (B), PFI (C) and DFI (D). E–G KM curves of CYP2C8 differentiating OS (E), DFI (F) and PFI (G) in selected cancer types

Fig. 13

Pan-cancer functional analysis and expression validation of CYP2C8 (A-D) Correlation of CYP2C8 expression with pan-cancer immune cell infiltration (A), immune checkpoint expression (B), immunoinflammation-related functions (C) and stemness score (D). E Immunohistochemical picture of CYP2C8 in some cancer types

Discussion

BLCA, a prevalent and highly aggressive malignancy, significantly undermines patients' overall well-being [43]. It is of paramount importance to develop effective prognostic tools that can guide treatment decisions and optimize patient outcomes in BLCA. Recent endeavors have introduced numerous prognostic frameworks for BLCA, encompassing those reliant on genes tied to propionate metabolism [31], glycolysis [32], as well as ferroptosis- [33] and dendritic cell-related genetic signatures [44]. These frameworks have deepened our understanding of the prognostic trajectory of BLCA, laying a solid foundation for clinical interventions. Although ADME genes are recognized for their prognostic relevance in renal cancer [24, 25], non-small-cell lung cancer [28] and head and neck squamous cell carcinoma [27], their prognostic utility in BLCA is yet to be confirmed, showing restricted applicability in this context. Our investigation introduces a novel paradigm that sheds new light on prognostic evaluation and clinical management approaches for BLCA patients.

Bioinformatics, as a sharp edge of scientific research in the information age, provides many facilities and techniques for researchers to explore in depth. In this regard, we are greatly inspired by the scientific content of Liu's lab. For example, discovering the role of single gene in a single tumor, they argued in detail the identity of SCN3B and CDK2 in glioma [36, 38], and proved the role of AIMP1 [40], CNIH4 [41] in head and neck squamous cell carcinoma. As for the regulation of gene sets in single cancers, Liu et al. collated the prognostic potential of CD8 T-cell markers in breast cancer [34], emphasizing the impact of the IGFBP family in gliomas [45]. In addition, single-gene pan-cancer studies have also inspired us, such as the discovery of the broad diagnostic and prognostic value of RAD51 [39], TRPM7 [46], and CENPA [47] in pan-cancer. Finally, the pan-cancer analysis of gene sets is more prospective, which will be the focus of our future research. In this regard, voltage-gated sodium channels [48] and disulfidptosis-related genes [49], have a wide range of regulatory effects in pan-cancer, which are worth exploring in depth.

Our analysis revealed that over a quarter (80/298, 26.8%) of ADME genes display aberrant expression patterns in tumor tissues, hinting at their crucial role in cancer initiation and progression. Employing univariate Cox regression and lasso regression methodologies [45], we pinpointed 17 ADME genes with prognostic implications: CYP2C8, SLC22A1, ABCA4, ABCC9, CHST6, CYP1B1, CYP27B1, CYP3A43, CYP3A7, CYP4F8, NR1I2, PON1, SLC22A11, SLC22A3, SULT1A3, TAP1, and UGT2B4. Current research underscores the involvement of ADME genes in tumorigenesis. Notably, CYP2C8, a pivotal element in our constructed PPI network, was conspicuously upregulated in BLCA tissues. Qu et al. notes a pronounced correlation between CYP2C8 polymorphisms and heightened bladder cancer vulnerability [50]. Moreover, CYP2C8, a key gene for iron metabolism, is closely related to the prognosis as well as the treatment of BLCA patients [51]. Additionally, transport proteins like SLC22A1, SLC22A3, and SLC22A11, belonging to the organic cation transporter (OCT) family, play a critical role in modulating the toxicity and effectiveness of platinum-based chemotherapy treatments [52]. ABC transporter superfamily members, ABCA4 and ABCC9, showing the highest mutation rate in the 17 prognostic genes, have been implicated in breast cancer prognosis and drug response [53, 54]. CHST6, a carbohydrate sulfotransferase-encoding gene, is associated with macular corneal dystrophy [55, 56] and may function as a prognostic indicator for low-grade gliomas [57]. NR1I2, a xenobiotic-sensing nuclear receptor, halts liver cancer progression by inhibiting EMT in hepatocellular carcinoma [58]. PON1 has emerged as a prognostic biomarker for hepatocellular carcinoma [59]. TAP1, critical for MHC class I antigen processing and display, exhibits elevated expression across multiple cancer types and holds significant prognostic value [60]. Variations in the UGT2B4 gene are linked to an increased likelihood of developing pancreatic cancer and influence the pathological outcomes of chemotherapy in patients with breast cancer [61, 62]. CYP1B1 promotes metastatic cancer cell growth in a fatty acid-dependent manner [63]. CYP27B1, vital for steroid biosynthesis, curbs ovarian cancer cell proliferation, migration, and invasion [64]. The CYP3A family contributes to the transition from chronic atrophic gastritis to gastric cancer and modulates chemotherapeutic metabolism and resistance [65]. Lastly, dopamine-mediated upregulation of SULT1A3 markedly enhances the metastatic potential of hepatocellular carcinoma cells [66].

We harnessed these 17 ADME genes associated with prognosis to devise a risk assessment metric for each specimen and constructed a predictive framework. In both the training and validation groups, individuals identified as high-risk and low-risk exhibited markedly different survival paths. Demonstrated by ROC curves, the efficacy of our model to predict overall survival was convincingly established. In addition, the model's accurate predictive effect was verified again in cohorts with different clinical characteristics. Furthermore, both univariable and multivariable Cox regression analyses independently confirmed the reliability of our risk metric as a predictor of outcomes in bladder cancer. Additionally, we crafted a nomogram that seamlessly integrates the ADME-based risk metric with clinical indices. And calibration, DCA and ROC curve analyses substantiate the exemplary prognostic performance showcased.

In addition, the choice of tools for assessing the prognostic ability of the model deserves an in-depth discussion. The Kaplan–Meier (KM) survival curve is a widely used nonparametric method for survival analysis, which can visualize differences in survival among different groups, but has some limitations. The first thing to note is that the KM curve is only suitable for univariate group comparisons, so when we need to include covariates in the analysis, Cox regression models can be utilized, which is also reflected in our study. Secondly, when the sample size included in the analysis is small, the KM curves are assessed to be unstable as the confidence intervals become wider. This situation can be considered by using Bayesian methods to introduce prior information to enhance the estimation stability. In addition, the KM curve analysis assumes that censoring is not related to event risk, but in reality, there may be informative censoring due to deterioration of the condition, which may result in biased results. We can try to perform a sensitivity analysis to assess the effect of the censoring mechanism on the results, or construct a competing risk model for assessment when there is competition for events (e.g., death or relapse). In addition to the KM Survival Curve, both the Receiver Operating Characteristic (ROC) Curve and Decision Curve Analysis (DCA) were implemented in this study to assess model capability, with different and complementary focuses.The ROC Curve complements each other by Area Under Curve ( AUC) to quantify the discriminatory power of the model, which is very intuitive and can also be combined with survival time to reflect the predictive power of the model at different time points. However, the AUC only reflects the statistical performance of the model and fails to take into account the actual needs of clinical decision-making, which is complemented by the DCA, which evaluates the clinical utility of the model by calculating the net benefit, and determines whether the model's use is superior to the default strategy of all-or-nothing intervention. At the same time, however, DCA ignores the accuracy of the model, so calibration curves or ROC curves are often co-opted to verify that the predicted probabilities are consistent with the actual risks. Co-using ROC curves and DCA realizes the unity of science and practicality and advances the clinical translation of the model.

Our study examined the functional annotations of genes that were expressed differently between risk groups, aiming to elucidate the role of ADME genes in BLCA. Gene Ontology enrichment analyses illuminated an overrepresentation of functions tied to the extracellular matrix (ECM) and signaling cascades, embracing ECM organization, extracellular structural architecture, collagen-enriched ECM frameworks, signaling receptor activation, and receptor-ligand interactions. The tumor microenvironment (TME) relies fundamentally on ECM, which is made up of a complex variety of proteins secreted by cells. These proteins are essential for maintaining structural integrity and enabling detailed cellular interactions [67]. Perturbations in ECM levels can directly alter its biological attributes, thereby fostering cancer dissemination via tumor cell heterogeneity [68]. KEGG demonstrates enrichment for many notorious pro-oncogenic pathways, such as PI3K-AKT signaling pathway and focal adhesion. Over-regulation of the PI3K-AKT pathway tends to mediate malignant cell proliferation and promotes the manifestation of cancer phenotypes [69]. In addition to mediating cell adhesion and promoting cancer cell proliferation, adhesion plaques also have significant effects on stemness regulation and cell proliferation [70]. GSEA, on the other hand, highlights the importance of immune-related pathways between high and low risk groups.

Our exhaustive investigation unveiled disparities in immune cell infiltration patterns across risk groups. Our observations highlighted a notable surge in M2 macrophages within the high-risk cohort, with the risk metric directly mirroring the abundance of all macrophage subtypes. Converging evidence underscores the association between elevated tumor-associated macrophages (TAMs) and unfavorable survival outcomes, as well as compromised therapeutic responses to chemotherapy and BCG therapy [71,72,73]. Notably, M2 macrophages are implicated in dampening immune responses and fostering angiogenesis in bladder cancer, thereby expediting tumor progression and dissemination [74]. This analysis indicates that the presence of M2 macrophages correlates with poorer outcomes among individuals in the high-risk category.In the group identified as high-risk, a pronounced deficit was observed in both activated dendritic cells and natural killer cells. Notably, the dendritic cell count inversely correlated with the risk metrics, underscoring a substantial imbalance in immune cell distribution. Given the pivotal role of dendritic cells in activating NKT and γδT cells to combat bladder cancer cells [75], the diminished occurrence of these cells may be a factor in the poorer outcomes observed in the high-risk cohort. Contrary to this, elevated levels of CD8 T cells in the low-risk group align with previous findings [76], suggesting that greater infiltration of these cells improves patient outcomes. And when we focused on the regulation of immune infiltration by individual model genes, it has been reported that the peroxidase-like activity of PON1 inhibits the differentiation of monocytes to macrophages [77]. In our model, PON1 expression was higher in the low-risk group, and cibersort analysis showed less macrophage infiltration in the low-risk group, and PON1 expression was negatively correlated with M0 macrophages, in agreement with the above literature. Moreover, TAP1 expression promotes cytotoxic T cell infiltration [78], whereas TAP1-deficient mice exhibit reduced expression of HLA class I molecules and thus a severe lack of CD8 T cells [79]. This is consistent with the higher expression of TAP1 as well as higher abundance of CD8 T cells in the low-risk group in our model. In addition to the above, there are fewer studies on the mechanism of model genes and immune cell infiltration. Although our study revealed a strong association between ADME genes and immune cell infiltration, the specific regulatory mechanisms remain unclear. Regulation of individual model genes Given the complexity of the immune microenvironment, the effect of individual genes may not fully explain the extent of immune cell infiltration, so in the future, we will try to conduct more precise controlled variable experiments and obtain samples from high- and low-risk patients from the clinic for detailed mechanism analysis.

Exploring correlations uncovered a pronounced inverse association between the risk metric and several immune checkpoint genes, such as CD40, TNFRSF14, TNFRSF25. CD40, a familiar immune checkpoint, is abundantly expressed on antigen-presenting cells, macrophages, and even tumor cells, where it interacts with helper T cells and mediates the enhancement of anti-tumor T cell responses [80]. In addition, it has been reported that CD40 stimulation of macrophages directly mediates cancer cell killing as well as microenvironment collapse [81]. The negative correlation index of TNFRSF14 and TNFRSF25 with risk scores was as low as −0.4. TNFRSF14 was associated with a better prognosis through promotion of apoptosis in BLCA [82], whereas TNFRSF25 mediated toxic killing of CD8 T cells to exert antitumor immunity [83]. Combined with our findings of upregulation of classical HLA expression and functional upregulation of MHC class I in patients in the low-risk group, the feasibility of diversified immunotherapies is highlighted. Our subsequent research illuminated a more favorable response to immunotherapy among low-risk patients within the IMvigor210 cohort. Our study revealed marked disparities in the tumor immune microenvironment among different risk categories, highlighting its importance for prognostic evaluations and the clinical management of bladder cancer. While focusing on specific model component genes, we found that fewer studies on the exact mechanism were reported. Among them, TAP1 was identified as a reliable prognostic target for pan-cancer immunotherapeutic response [60], which increases PDL1 expression and can be enhanced by oxaliplatin stimulation [84]. In addition, in colorectal cancer, Chen et al. found that CYP1B1 could trigger anti-PD-1 immunotherapy inactivation by mediating ubiquitination degradation of ACSL4 [85]. In conclusion many ADME genes may play a potentially influential role in cancer immunotherapy. Future endeavors should delve deeper into the intricate interplay among ADME genes, immune cell dynamics, and the immunotherapeutic response, with the aim of enhancing patient outcomes.

Malignant phenotypic correlation analysis showed that our established risk score showed a highly positive correlation with angiogenesis as well as EMT, which contributes to an in-depth understanding of the root causes of the progression of poorer tumor prognosis. Interestingly, we found that TMB showed a significant negative correlation with the risk score, implying that the low-risk group had a higher propensity for mutation. Previous study [86] have elaborated that antigenic changes due to mutations promote immune cell recognition and presentation, which mediates cell killing, resulting in greater benefit from immunotherapy. This is another aspect that confirms the better prognosis of immunotherapy in low-risk patients that we found in the above study.

To refine clinical outcomes, we performed a drug sensitivity assessment, revealing distinct response patterns between risk groups. Our research reveals that patients with high-risk profiles generally exhibit positive responses to traditional medications such as Foretinib and Dasatinib. Conversely, these same groups show resistance to drugs like Oxaliplatin, Gefitinib, and Cyclophosphamide. When we looked at the mechanism by which specific ADME model genes affect chemotherapy response, it had been reported that overexpression of CYP1B1 mediates complex cancer cell phenotypes such as enhanced mesenchymalization, which led to resistance to chemotherapeutic agents such as 5-fluorouracil, cisplatin, and others [87]. In addition, platinum-based chemotherapeutic agents such as oxaliplatin had been reported to be substrates for SLC22A1-3 [88], and high expression of SLC22A3 promoted the uptake of oxaliplatin, which accumulated more platinum and exerted stronger toxic effects on cancer cells [89]. However, in our study, the trend of differential expression of SLC22A1-3 in the high- and low-risk groups varied, so the specific regulatory mechanism threw further control variables need to be further controlled variables in-depth study. Moreover, the use of chemotherapeutic agents also had an impact on the activity of the ADME gene protein, such as Gefetinib and other kinase inhibitors, which were inhibitors of the activity of CYP3A [90]. We speculated that the reduction of the activity of the ADME protein may also subsequently regulate the efficacy of the chemotherapeutic agents in a reverse manner. There are few reports on the regulatory mechanisms of specific ADME genes in chemotherapy, so future research should focus on this aspect to further advance customized treatment strategies and combination chemotherapy regimens.

CYP2C8 was extracted as a core target by us, and by qRT-PCR with immunohistochemistry, we perfectly validated the differential expression conclusion obtained from bioinformatics analysis. In our analysis, high expression of CYP2C8 predicted a better prognosis for BLCA patients, suggesting that CYP2C8 may mediate some protective mechanisms. Notably, high expression of CYP2C8 was associated with higher infiltration of dendritic cells, which, as the mainstay of antigen presentation, have been reported to positively contribute to tumor immunity in bladder cancer [76]. We speculate that CYP2C8 may contribute in regulating dendritic cells, mediating a better prognosis by modulating the immune microenvironment. GSVA heatmap showed that the high expression group exhibited higher peroxisome activity. One of the major functions of peroxisomes is to metabolize reactive oxygen species [91], which act as active second messengers in the cell and are involved in proliferation, angiogenesis, and other notorious pro-carcinogenic processes [92]. We speculate that CYP2C8 may mediate peroxisome-associated functions to stabilize the intracellular reactive oxygen species environment, thereby reducing the malignancy of cancer cells to a certain extent and leading to a better prognosis. In addition, there were significant metabolic pathway differences between the high and low expression groups of CYP2C8, which may also have an impact on cancer development. We then saw the light and pushed CYP2C8 from bladder cancer to pan-cancer and revealed its potential as a novel prognostic target through comprehensive immune-related and stemness-related analyses. Although CYP2C8 is closely associated with the prognosis of bladder cancer and even pan-cancer, the specific regulatory mechanism remains unclear. In the future, we will focus on the identity of CYP2C8 in the metabolic microenvironment and immune microenrivronment of cancer cells and explore the specific mechanism by which it affects cancer progression.

Personalized medicine is an increasingly popular model of care for clinicians, aiming to provide the most rational healthcare plan by integrating the patient's multi-omic data, clinical characteristics, physiological conditions, living environment and even economic status [93]. Personalized medicine is particularly important in the management of BLCA, a urologic tumor with significant risk heterogeneity, where cancer progression, genomic landscapes, and immune microenvironments vary from patient to patient [94, 95]. With an increasing number of drug candidates for BLCA [96], the selection of rational therapies is a pressing issue [97], and by utilizing the differences in drug sensitivities obtained from the risk stratification of our model, clinicians can attempt to formulate precise and targeted therapeutic strategies. For patients identified as high-risk, they may show higher sensitivity to Staurosporine, Dasatinib, and thus these drugs can be selected as first-line treatment options. In addition, for the treatment of patients in the high-risk group, it is necessary for clinicians to work together in a joint multidisciplinary team (MDT) because the comprehensive management of patients with BLCA, including physical, nutritional, and emotional aspects, is considered to be a crucial part of the treatment process [98]. On the other hand, high-risk patients may show a tendency to have a worse prognosis, which highlights the importance of choosing the desired outcome [99] and may require the cooperation of the medical staff for appropriate end-of-life care. For patients identified as low risk, according to our model, may show higher sensitivity to anti-PDL1 immunotherapy as well as Oxaliplatin, 5-Fluorouracil chemotherapy. Risk modeling should be used for ongoing prognostic assessment while clinicians administer targeted therapies based on protocols. When periodic pathology and tumor marker assessments determine that the tumor is likely to progress, clinicians can use risk modeling to rapidly reassess the patient's changing risk of regression and adjust treatment strategies accordingly. With our risk model, clinicians have a reference point to implement more precise treatments rather than relying on empirical or non-specific approaches, thus facilitating the process of personalized management of BLCA and laying a solid foundation for ensuring a better prognosis for patients.

While our findings are encouraging, several limitations must be acknowledged. The reliance on publicly available datasets (mainly TCGA and GEO) may introduce biases related to sample selection and population characterization, which will lead to a lack of model generalizability within large groups. For example, the TCGA database used for model training focuses on European and American populations and lacks ethnic diversity such as Asian and African populations, leading to biases in genetic background, environmental exposure, and other factors. Moreover, the samples are mainly early-stage cancer patients, with relatively few advanced and metastatic cases, which limits the applicability of the results to advanced cancers to a certain extent. And the GEO dataset is mainly expression profiles from high-throughput sequencing, which is not a comprehensive record of patients' clinical characteristics, drug treatment response, immunotherapy effect and other parameters, resulting in a lack of clinical efficacy exploration of the model. Therefore, we are trying to construct a more comprehensive BLCA cohort, including the inclusion of patients' treatment information, to address these limitations, further explore the generalizability of the model, and lay the foundation for more detailed mechanistic analysis. In addition, technical and biological biases in bioinformatics analyses are important factors contributing to the limitations [100, 101]. Limitations in technical factors are mainly caused by sequencing methods, such as for microarray chips microarray chips, where similar sequences may hybridize with the wrong probes, leading to false signal detection or signal interference. And RNA sequencing may have problems such as PCR amplification bias and differences in sequencing depth and coverage. To cope with these limitations, we can take a specific approach to specific data, such as analyzing microarray data with the “limma” package, and analyzing RNA sequencing data with “DESeq2” or “edgeR” packages. RNA sequencing data are analyzed with “DESeq2” or “edgeR”, and the sequencing depth, batch effect and amplification bias are adjusted by normalization procedures. It is also worth noting that cross-platform batch effects can mask the true biological signal, so when combining datasets from multiple platforms for analysis, the “sva” package can be used for de-batching. Limitations due to biological factors are mainly due to intra-tumor heterogeneity, sample purity, etc. Bulk RNA-seq averages the signals from all cells in the sample and is not able to show the intra-tumor expression landscape at a fine resolution. Therefore, the detailed molecular landscape can be further explored in the future by means of single-cell sequencing and other means. Finally, the validation of the experiments is a very direct guide to the bioinformatics findings, and more in-depth mechanistic experiments should be conducted in the future to explore the role of ADME genes in the progression of BLCA.

Conclusion

Our research established and confirmed a prognostic risk assessment tool based on ADME gene utilization, which showed high predictive accuracy for survival rates in bladder cancer. Functional enrichment analyses revealed unique immune-related signatures across high- and low-risk patient groups, connecting ADME gene expression patterns with the complex dynamics of the tumor microenvironment. Furthermore, the integration of clinical factors into a nomogram serves as a valuable instrument for precise prognosis prediction, enabling tailored treatment strategies for BLCA patients, thereby advancing personalized medicine. In addition, we identified CYP2C8 as a novel BLCA target and generalized it to pan-cancer.