Tumor Immune Infiltration as a Predictor of Response to Neoadjuvant Chemo-Immunotherapy in Muscle-Invasive Bladder Cancer: An Integrative TCGA Analysis

1. Introduction

The muscle-invasive bladder cancer (MIBC) is still among the most aggressive types of urologic malignancies, with high rates of recurrence, early metastasis, and high mortality even with improvements in the treatment [1]. The conventional management entails neoadjuvant chemotherapy (NAC) containing cisplatin followed by radical cystectomy, which offers survival advantages but is associated with substantial variability in pathological response, including inconsistent rates of tumor downstaging and pathologic complete response [2,3]. Over the last several years, immune checkpoint inhibitors have transformed the management of advanced bladder cancer, leading to growing interest in perioperative strategies, including approaches in which immunotherapy is administered alongside neoadjuvant chemotherapy and, in selected studies, continued in the adjuvant setting [4]. The combination of the two strategies is intended to increase tumor antigen release, immune priming and induce long-lasting antitumor immune responses [5].Nevertheless, responses to chemo-immunotherapy remain heterogeneous across clinical trials; for example, the NIAGARA trial reported that the addition of the PD-L1 inhibitor durvalumab to neoadjuvant chemotherapy improved pathologic complete response rates by approximately 10%, highlighting that only a subset of patients derive substantial benefit [6]. Determining predictive biomarkers that can be used to differentiate between responders and non-responders prior to the commencement of treatment has thus become a major concern [7].

The Cancer Genome Atlas (TCGA) has produced both multi-omics datasets of bladder cancer, which has provided an unparalleled chance to study molecular properties that could be used to regulate the response to therapy [8]. MIBC has been reported to have high genomic diversity, including the changes in DNA damage response (DDR) genes, oncogenic mutations, and dynamic tumor mutation burden (TMB) [9,10]. Such properties dictate tumor behavior, immune response and sensitivity to cytotoxic therapy, implying that they might prove useful as biomarkers of response to perioperative chemo-immunotherapy [11]. For example, alterations in DNA damage repair (DDR) pathways have been associated with increased sensitivity to platinum-based chemotherapy and may also augment genomic instability and neoantigen load, potentially enhancing tumor immunogenicity and responsiveness to immune checkpoint blockade [12]. On the same note, elevated TMB and distinct signature of gene expression related to interferon signaling or infiltration by cytotoxic lymphocytes have demonstrated potential as predictors of immunotherapy efficiency on other types of malignancies [13].

Besides genomic changes, the tumor microenvironment is the primary element in the regulation of the response to treatment [14]. The patterns of immune infiltration, including the presence of CD8+ T cells, activated dendritic cells, and tumor-associated macrophage may be used to identify a tumor as an inflamed, immune-excluded, or immune-desert [15]. These differences directly relate to the effectiveness of immune-based therapy. There is therefore a need to combine genomic, transcriptomic and immune-cell profiling in developing a holistic approach to understanding the biological determinants of treatment outcome [16,17]. Nevertheless, there are limited data that directly associate such biomarkers with a response to perioperative chemo-immunotherapy in MIBC [18].

Due to the accelerated rate of combination neoadjuvant trials, predictive biomarkers are in high demand to inform individual-based treatment and prevent unwarranted toxicity imposed by futile therapies [19]. The use of TCGA-BLCA data enables scientists to research the molecular landscape of MIBC on a large scale, define potential biomarkers based on the real world of genomic heterogeneity, and establish a mechanistic model of predicting therapeutic response [20]. The objective of this study is to identify genomic and immune-related biomarkers, including gene expression profiles, mutation and DNA damage response status, tumor mutational burden, and immune infiltration signatures, within the TCGA-BLCA cohort that characterize immune-active tumor phenotypes and may be associated with responsiveness to neoadjuvant chemo-immunotherapy in muscle-invasive bladder cancer.

The perioperative treatment landscape for muscle-invasive bladder cancer is evolving beyond surgery alone, with increasing integration of systemic therapy earlier in the disease course [1,4,18,19]. Although cisplatin-based neoadjuvant chemotherapy remains standard and achieves complete responses in a subset of patients [2,3], chemotherapy can also enhance antitumor immunity through immunogenic cell death, supporting the rationale for combination strategies with immune checkpoint inhibitors [5,6]. Ongoing perioperative trials further reflect the shift toward earlier incorporation of immunotherapy in curative-intent settings [4,6,18,19]. At the same time, substantial effort has focused on identifying predictive biomarkers, including tumor mutational burden, DNA damage response alterations, immune infiltration signatures, and molecular subtypes characterized through TCGA and related genomic platforms [8,9,10,11,12,13,14,15,16,17,20,24]. These developments emphasize the need for integrative genomic and immune-based analyses to better characterize tumor phenotypes that may inform selection for perioperative chemo-immunotherapy, providing a clear rationale for the present TCGA-based investigation.

These insights can be used to hasten precision oncology in bladder cancer to inform clinical trial stratification and can aid in improving personalized decision-making in patients with MIBC.

2. Methodology

2.1. Study Design and Data Source

This study was designed as a retrospective, bioinformatics-based analysis utilizing publicly available data from The Cancer Genome Atlas Bladder Urothelial Carcinoma (TCGA-BLCA) project [21]. TCGA provides harmonized genomic, transcriptomic, and clinical datasets derived from treatment-naive primary tumor specimens. Data extraction and preprocessing were performed using R statistical software (version 4.5.0) with the TCGA biolinks package, which enables standardized retrieval and integration of TCGA molecular and clinical data. Because TCGA does not include direct information on response to neoadjuvant chemo-immunotherapy, this analysis focused on biologically and clinically relevant surrogate biomarkers that have been previously associated with immunotherapy sensitivity in bladder cancer and other solid tumors. These included immune infiltration metrics derived from transcriptomic data, tumor mutational burden, and alterations in DNA damage response genes, which have been linked to tumor immunogenicity, neoantigen load, and likelihood of benefit from immune checkpoint blockade.

2.2. Study Population

The study population consisted of patients with histologically confirmed bladder urothelial carcinoma included in the TCGA-BLCA cohort. Samples were eligible for inclusion if complete RNA sequencing data, somatic mutation data, and key clinical variables (age, sex, and pathologic stage) were available. To ensure analytic consistency, samples with missing data for key variables required for immune profiling, genomic characterization, or multivariable modeling were excluded. All analyses were conducted at the patient level using harmonized TCGA identifiers to prevent duplication and misclassification.

2.3. Variables and Measure

Gene expression profiling was conducted using RNA sequencing data normalized as fragments per kilobase of transcript per million mapped reads, upper quartile (FPKM-UQ). Expression data were processed to generate gene-level matrices by removing Ensembl version numbers and mapping Ensembl identifiers to official gene symbols. Only uniquely mapped genes were retained, and duplicated symbols were excluded to prevent redundancy. Expression values were log2-transformed to improve variance stability and cross-sample comparability.

Tumor immune infiltration was estimated using the xCell algorithm, a gene signature–based method that infers relative enrichment of immune and stromal cell populations from bulk RNA sequencing data. Immune variables of primary interest included ImmuneScore, CD8⁺ T-cell enrichment, and natural killer (NK) cell enrichment, selected for their biological relevance to antitumor immunity.

Somatic mutation data were derived from TCGA mutation annotation format files. Alterations in DNA damage response (DDR) genes were identified using a predefined DDR gene panel, and tumors were classified as DDR-mutated or DDR-wildtype. Tumor mutational burden was calculated as nonsynonymous mutations per megabase and dichotomized into high and low groups using the cohort median. Clinical variables included age, sex, and pathologic stage, which was grouped into stages I–II and III–IV for use as a covariate in multivariable analyses to account for differences in tumor biology. This categorization was not intended to imply eligibility for neoadjuvant treatment but rather to adjust for baseline disease severity within the TCGA cohort. Due to the lack of direct treatment response data, a binary surrogate outcome was constructed using immune-related transcriptomic features and genomic characteristics. This surrogate was designed to primarily reflect predicted likelihood of response to immunotherapy and chemo-immunotherapy combinations, rather than sensitivity to chemotherapy alone. This approach is consistent with established practices in translational immuno-oncology research and enabled stratification of tumors into “likely responder” and “non-likely responder” groups for comparative and predictive modeling analyses.

2.4. Statistical Analysis

Baseline demographic, clinical, genomic, and immune characteristics were summarized by inferred response status using descriptive statistics. Continuous variables were reported as mean (standard deviation) and compared using Student’s t-tests, while categorical variables were presented as counts (percentages) and compared using Pearson’s chi-square tests. Multivariable analysis was performed using Firth’s penalized logistic regression to identify independent predictors of response, accounting for separation and small-sample bias. Predictor variables included immune infiltration metrics (ImmuneScore, CD8+ T cells, and NK cells), DDR mutation status, TMB group, age, sex, and stage group. Continuous predictors in the model were standardized using z-scores. Odds ratios with 95% confidence intervals were reported. Multicollinearity was assessed using variance inflation factors, which indicated no collinearity concerns. All analyses were conducted in R version 4.5.0.

2.5. Missing Data

Only samples with complete data for all variables required in descriptive and multivariable analyses were included. Missing data were not imputed, and complete-case analysis was applied to preserve the validity of statistical inferences. The extent of missingness was assessed prior to analysis to ensure adequate sample size and representativeness.

2.6. Ethical Considerations

All data used in this study were obtained from TCGA, a publicly available, de-identified resource. As no protected health information was accessed and no human subjects were directly involved, institutional review board approval and informed consent were not required in accordance with applicable guidelines.

3. Results

Table 1 presents the baseline demographic, genomic, and immune characteristics of patients from the TCGA-BLCA cohort stratified by inferred immune response status (non-likely responders versus likely responders).

From the baseline characteristics summarized in Table 1, the mean age of non-likely responders (67.80 ± 11.06 years) was comparable to that of likely responders (68.15 ± 9.68 years), with no statistically significant difference observed (p = 0.707), indicating that age was unlikely to confound immune response stratification. In contrast, sex distribution differed significantly between groups, with females representing a greater proportion of likely responders compared with non-likely responders 91 (35.0%) vs. 56 (22.3%), suggesting a potential sex-related influence on immune responsiveness. Pathologic tumor stage was similarly distributed, with stage III–IV disease present in 113 (43.5%) likely responders and 98 (39.0%) non-likely responders, and no significant difference detected (p = 0.355). DDR mutation status also did not differ significantly, as DDR-mutated tumors comprised 97 (37.3%) of likely responders and 105 (41.8%) of non-likely responders, indicating that DDR alterations alone did not distinguish immune response categories at baseline.

Tumor mutational burden demonstrated a non-significant trend toward enrichment among likely responders (p = 0.073), with a higher proportion of TMB-high tumors observed in this group ,116 (50.4%) vs. 98 (41.7%). In contrast, immune-related features showed pronounced differences. Likely responders exhibited higher ImmuneScore (0.23 ± 0.18 vs. 0.04 ± 0.02) and StromaScore (0.08 ± 0.10 vs. 0.05 ± 0.08), along with greater enrichment of CD8⁺ T cells, natural killer cells, and macrophages. Overall, these findings emphasizes the strong association between immune infiltration patterns and inferred responsiveness to neoadjuvant chemo-immunotherapy, as further illustrated in Figure 1.

Tumors characterized by both DDR mutations and high tumor mutational burden(TMB) demonstrated the highest ImmuneScore values, indicating a more immune-infiltrated tumor microenvironment. In contrast, tumors with low TMB, regardless of DDR status, exhibited lower immune infiltration. These findings reinforce the association between genomic instability and immune activation and provide biological context for the elevated immune metrics observed among likely responders in Table 1.

Table 2 shows the multivariable Firth’s penalized logistic regression model assessing the independent variables of likely immune response to neoadjuvant chemo-immunotherapy in muscle-invasive bladder cancer(MIBC).The model comprises of immune infiltration factors, genomic features and clinical covariates.

After multivariable adjustment, ImmuneScore emerged as the only independent predictor significantly associated with likely immune response(β = 131.00, p = 0.003).In contrast, individual immune cell populations, including CD8⁺ T cells (β = 0.42, p = 0.771) and natural killer cells (β = 0.23, p = 0.809), were not independently associated with response after accounting for overall immune context, suggesting that composite immune activity may be more informative than isolated cell subsets.

Genomic features, including DDR mutation status and tumor mutational burden, did not demonstrate significant independent associations with immune response in the adjusted model. Similarly, traditional clinical variables such as age, sex, and pathologic tumor stage were not significant predictors. Overall, these findings indicate that global immune infiltration, rather than individual immune cell types, genomic alterations, or clinicopathologic factors alone, independently drives inferred responsiveness to neoadjuvant chemo-immunotherapy in this cohort. To facilitate clinical interpretation of these regression coefficients, odd ratios with 95% confidence interval were calculated as shown in Table 3.

From the odds ratio findings its only ImmunoScore that demonstrated a statistically significant association with likely immune response(p=0.003),a clear indication of a strong positive effect with increased odds of response as overall immune infiltration increased. This study findings illustrates that global tumor immune infiltration as the core independent predictor of inferred responsiveness to neoadjuvant chemo-immunotherapy in this analytic group.

Collectively, these results suggests that even if multiple immune and genomic factors may lead to tumor ,overall immune activation within the tumor microenvironment is the primary driver of immune responsiveness in muscle-invasive bladder cancer.

4. Discussion

This study provides insight into the biological features associated with inferred responsiveness to neoadjuvant chemo-immunotherapy in muscle-invasive bladder cancer (MIBC). Consistent with an expanding body of evidence, the findings suggest that treatment responsiveness in MIBC is primarily driven by the tumor immune microenvironment rather than isolated genomic alterations. Across multiple analytic layers, global immune infiltration emerged as the most robust and independent determinant of predicted immune responsiveness.

An important context for interpreting these results is the well-established molecular heterogeneity of MIBC. Intrinsic molecular subtypes, particularly basal/squamous and luminal phenotypes, exhibit distinct immune and clinical behaviors [8,9,24]. Basal/squamous tumors are characterized by higher immune infiltration, increased PD-L1 expression, and enrichment of interferon-related signaling pathways features that align closely with the high ImmuneScore phenotype observed in this study. In contrast, luminal papillary tumors typically display lower immune infiltration and are often driven by FGFR3-associated signaling pathways, which have been linked to relative resistance to immunotherapy [24]. These subtype-specific immune landscapes provide a biologically plausible explanation for why global immune activation, rather than individual genomic alterations, dominated as the strongest predictor of inferred response in this study analysis.

In baseline comparisons, tumors classified as likely responders demonstrated significantly higher ImmuneScore, accompanied by increased infiltration of CD8⁺ T cells, natural killer (NK) cells, macrophages, and stromal components. These observations are concordant with prior studies demonstrating that an inflamed tumor microenvironment is closely associated with improved outcomes following immunotherapy-based regimens in bladder cancer and other solid malignancies [6,14,15]. The enrichment of cytotoxic and innate immune populations suggests the presence of an active antitumor immune milieu that may be further potentiated by chemotherapy-induced immunogenic cell death, enhanced antigen presentation, and amplification of immune checkpoint blockade efficacy [5].

Although PD-L1 expression and interferon-gamma (IFN-γ) related gene signatures were not directly assessed in this TCGA-based analysis, ImmuneScore likely captures these established biomarkers indirectly. Prior studies have demonstrated strong correlations between immune infiltration, PD-L1 expression, and IFN-γ driven chemokine signatures such as CXCL9 and CXCL10 [7,16,22]. Notably, neoadjuvant immunotherapy trials, including PURE-01 and ABACUS, have shown that pre-existing immune activation and PD-L1 positivity are among the strongest predictors of pathologic complete response [6,18,19,23]. The strong performance of ImmuneScore in this study therefore supports its use as a biologically meaningful surrogate marker reflecting convergent immune activation pathways known to mediate response to immune checkpoint inhibition.

While DNA damage response (DDR) gene alterations and higher tumor mutational burden (TMB) showed trends toward enrichment among likely responders, neither variable retained independent significance in multivariable analysis. This finding is noteworthy given prior reports linking DDR deficiency and elevated TMB to increased neoantigen generation and immunotherapy sensitivity [12,13].This study’s results suggest that these genomic features may contribute to immune activation indirectly, primarily through downstream modulation of the tumor immune microenvironment, rather than functioning as independent predictors of response. This interpretation is consistent with integrative multi-omics studies demonstrating that immune contexture frequently outweighs individual genomic alterations in predicting therapeutic outcomes [8,11].

Mechanistically, ImmuneScore likely outperforms single immune-cell metrics because it reflects the coordinated immune contexture of the tumor microenvironment. By integrating multiple dimensions of immune activity including antigen presentation, chemokine signaling, stromal–immune interactions, and leukocyte recruitment ImmuneScore captures a systems-level measure of immune readiness [15]. This may explain why individual immune cell populations such as CD8⁺ T cells and NK cells lost statistical significance when modeled alongside ImmuneScore, highlighting the limitations of relying on isolated immune parameters.

In the adjusted Firth’s penalized logistic regression model, ImmuneScore emerged as the sole independent predictor of inferred immune responsiveness, maintaining a strong positive association after accounting for genomic variables, immune subsets, and clinical covariates. These findings align with pan-cancer immune landscape analyses emphasizing the primacy of global immune activation over individual cellular or genomic features when predicting response to immunotherapy [15].

Collectively, these results support a biologically plausible framework in which overall immune activation within the tumor microenvironment is the principal driver of responsiveness to combined chemotherapy and immunotherapy in MIBC. This framework is summarized in the conceptual model shown in Figure 2.

This concept aligns with emerging neoadjuvant clinical trial data demonstrating that baseline immune infiltration and interferon-related gene expression signatures are more reliable indicators of pathologic response than genomic alterations alone [6,16,18]. As neoadjuvant chemo-immunotherapy becomes increasingly incorporated into clinical practice, these findings underscore the potential value of immune profiling in patient stratification and clinical trial design [4,19].

An additional observation of interest was the higher likelihood of inferred responsiveness among female patients. Emerging evidence suggests that sex-based differences in immune regulation may influence tumor–immune interactions in bladder cancer, potentially through hormonal effects and X-linked immune regulatory pathways [7]. While this finding should be interpreted cautiously given the retrospective nature of the dataset, it highlights an area warranting further investigation in prospective immunotherapy studies.

Clinical Implications, Strengths, Limitations of the Study

This study has important clinical implications. Immune-based biomarkers such as ImmuneScore may aid in identifying patients more likely to benefit from the immunotherapy component of neoadjuvant chemo-immunotherapy regimens, potentially helping to spare patients with immune-cold tumors from unnecessary exposure to checkpoint inhibitors and their associated toxicities [7,11]. Furthermore, high immune infiltration may help identify candidates for emerging bladder-sparing or organ-preservation strategies, particularly when integrated with radiologic assessment and circulating tumor DNA (ctDNA) monitoring [2,18].

This study utilized a large, publicly available, and well-curated dataset from TCGA-BLCA, enabling a comprehensive integrative analysis of genomic, immune, and clinical features in muscle-invasive bladder cancer. The use of standardized RNA-sequencing and mutation data ensures high data quality and reproducibility. A key strength is the incorporation of multiple immune infiltration metrics alongside DDR mutation status and tumor mutational burden, allowing for a biologically informed assessment of immune responsiveness. Additionally, the application of Firth’s penalized logistic regression strengthens inference in the presence of rare events and quasi-complete separation, improving the stability and reliability of effect estimates.

Several limitations should be acknowledged. First, TCGA lacks direct clinical information on neoadjuvant chemo-immunotherapy administration and treatment response; therefore, immune responsiveness was inferred using surrogate immune and genomic markers rather than observed clinical outcomes. These surrogate measures primarily reflect inferred immunotherapy sensitivity and do not directly capture chemotherapy benefit. In addition, chemotherapy itself can alter the tumor microenvironment through mechanisms such as enhanced antigen release, modulation of immune cell infiltration, and changes in stromal composition. Because these dynamic treatment-induced effects are not captured in baseline TCGA tumor samples, the interaction between chemotherapy-mediated microenvironment remodeling and immune-related biomarkers could not be directly assessed. Second, immune cell estimates were derived from bulk RNA-seq data rather than single-cell or spatial profiling, which may mask intratumoral heterogeneity. Third, multiple ongoing phase III perioperative trials (e.g., KEYNOTE-905, EV-304, VOLGA) are expected to provide definitive clinical evidence regarding the benefit of chemo-immunotherapy combinations and the role of biomarkers such as immune infiltration. Finally, some regression estimates particularly for ImmuneScore were influenced by quasi-complete separation, resulting in extremely large effect sizes that, though its clinically meaningful it should be interpreted cautiously.

Future research should validate these findings in independent cohorts with documented neoadjuvant chemo-immunotherapy exposure and clinical response data. Prospective studies integrating longitudinal sampling would help clarify how immune infiltration evolves with treatment. Incorporation of single-cell approaches could provide finer resolution of immune cell interactions within the tumor microenvironment. Finally, combining TCGA-based biomarkers with real-world clinical trial datasets may enhance translational relevance and support the development of clinically actionable immune-genomic prediction models.

5. Conclusions

This study evaluated genomic and immune-related predictors of inferred responsiveness to neoadjuvant chemo-immunotherapy in muscle-invasive bladder cancer using TCGA data. Global tumor immune infiltration emerged as the strongest predictor of presumed treatment response, whereas individual genomic alterations, including DNA damage response status and tumor mutational burden, were not independently predictive. Clinical variables did not significantly contribute to response stratification. These findings emphasize the importance of immune contexture over isolated genomic features in determining treatment sensitivity. Incorporating immune profiling may enhance patient stratification and inform the design of future neoadjuvant trials. Prospective studies incorporating direct treatment response data are needed to validate these observations and support their clinical application.