To date, immune checkpoint inhibitors (ICIs) have emerged as a leading treatment for metastatic cancer, significantly improving patient survival while causing relatively few side effects. However, the objective response rate for ICIs remains low approximately 30% in urothelial carcinoma (UC), underscoring the urgent need for predictive response biomarkers. Several state-of-the-art signatures have been revealed in top-tier journals, highlighting the importance of this field. As the number of genes (~20,000) far exceeds the sample sizes of typical training sets (generally ≤ 300), we first developed feature selection procedures to reduce the number of features to a few hundred. We then trained multiple machine learning classifiers using the selected genes and the IMvigor210 dataset, which includes RNA-seq and clinical data from ~298 patients with metastatic UC (mUC). Notably, our predictor LogitDA, using the identified 49-gene signature, achieved a prediction AUC of 0.75 in an independent dataset, PCD4989g(mUC). Moreover, our signature outperformed six state-of-the-art signatures, PD-L1 IHC, and five tumor microenvironment signatures, including IFN-γ, T-effector, and T-cell exhaustion signatures. When we integrated each of the six known signatures with our own, our signature still surpassed the integrated ones in terms of prediction AUC and accuracy in the PCD4989g(mUC) dataset. From our signature, we identified key prognostic biomarkers, with the top five markers LYRM1, RFC4, CENPL, SPAG5, and CACYBP (Benjamini-Hochberg adjusted P < 0.0025) in the IMvigor210 dataset. Finally, we performed pathway analyses using Reactome (MSigDB) and KEGG, to reveal some immune-related pathways enriched such as MHC class II antigen presentation.
biomarkercancerimmunotherapymachine learningregressionpredictionThe author(s) declare financial support was received for the research and/or publication of this article. This research was supported in part by Academia Sinica, Taiwan (Group research grant to GS); National Defense Medical Bureau (MND1214075 to Y-TT); Tri-Service General Hospital Research Foundation (TSGH_E_114231 to T-LC) and National Science & Technology Council, Taiwan, Republic of China (113-2118-M-001–002 and 114-2118-M-001-001 to GS to T-LC).section-at-acceptanceCancer Immunity and ImmunotherapyIntroduction
Metastasis accounts for nearly 90% of cancer-related deaths and remains a major challenge in effective cancer treatment. Immune checkpoint inhibitors (ICIs) have improved outcomes in several metastatic cancers. In metastatic urothelial carcinoma (mUC) (1–3), the PD-L1 inhibitor atezolizumab has shown durable clinical efficacy (4). Atezolizumab, a humanized monoclonal antibody, binds PD-L1 and blocks its interaction with PD-1 and B7.1, thereby restoring tumor-specific T-cell immunity. However, only about 20% of mUC patients achieve objective responses (5, 6), highlighting the need for reliable biomarkers to predict treatment benefit before therapy initiation (7).
Banchereau et al. reported that tumor mutation burden (TMB) and PD-L1 expression had limited predictive power across the IMvigor210 (mUC), POPLAR, and IMmotion150 (RCC) cohorts, whereas RNA-seq-based models captured their effects more accurately (6). They further demonstrated that cancer-specific models outperform pan-cancer approaches. Guided by this principle, we developed an mUC-specific transcriptomic signature predictive of response to atezolizumab, aiming to stratify patients most likely to benefit from ICI therapy.
Data-driven machine learning (ML) models often fail to generalize response predictions to independent datasets (8). To address this problem, we applied transfer learning to enhance feature selection and classifier training. Because the number of genes far exceeds the number of ICI-treated samples, effective feature reduction is critical. We incorporated unsupervised domain adaptation (DA) to reduce ~17,000 genes to a few hundred, prioritizing those with similar statistical distributions across training and test datasets (9–11). This ensured that selected features retained predictive power across domains.
Following feature selection, we performed cross-validation of four classifiers, logit, lasso, support vector machine (SVM), and random forest, using the IMvigor210 and IMmotion150 datasets to identify the most consistent predictor. As illustrated in Figure 1A, LogitDA, a logistic regression-based model, consistently achieved the highest prediction accuracy across both cohorts. Consequently, LogitDA and its derived gene signature were selected for independent validation in the PCD4989g(mUC) dataset.
The overall scheme, and predictive performance of our signature versus the six known ones in independent mUC dataset. (A) The flowchart of our scheme, where the DE genes of IMvigor210 were selected by FDR < 0.10. (B) The scheme of training and test of logit-based predictors with our signature and the six known ones. The datasets and signatures used to train and testing the logit model, and the number of samples in each dataset are shown. (C) The area under the receiver operating characteristic curve (AUC) for PCD(mUC) is displayed. The random expectation (AUC = 0.5) is shown in dotted lines.
Flowchart detailing the overall scheme and the predictive performance of our signature (Ours) versus the six known ones. Panel A outlines steps for identifying DE genes, and training the predictors. Panel B shows the training dataset (IMvigor210, n=298) and test dataset (PCD4989g mUC, n=94), highlighting logit-based predictors and features for training. Panel C displays ROC curves comparing two predictors with AUCs of 0.75 and 0.55, and a bar graph showing various AUCs for different signatures, with “Ours” scoring the highest.
We further compared our 49-gene mUC signature with six established immunotherapy response predictors, including PD-L1 IHC and five tumor microenvironment (TME)-associated gene signatures: the IFN-γ signature (1), T-cell dysfunction signature (12), T-effector (tGE8) signature (2), T-cell–inflamed GEP (13), TIDE T-exhaust signature (12), and CD8T signature (14). The IFN-γ and T-cell-inflamed GEP each comprise 18 genes linked to antigen presentation, chemokine signaling, and adaptive immune resistance, and have shown pan-cancer predictive value for anti-PD-(L)1 therapy. The tGE8 signature, composed of IFNG, CXCL9, CD8A, GZMA, GZMB, CXCL10, PRF1, and TBX21, correlates with PD-L1 expression and is elevated in responders (2, 15). The TIDE model identified 50 genes associated with T-cell exhaustion (12), while the CD8T signature (CD8A, CD8B, GZMA, GZMB, PRF1) reflects CD8+ T-cell activity (14).
We found that our 49-gene signature for mUC outperformed all six known signatures in terms of prediction AUC and accuracy in an independent test dataset PCD4989g(mUC). We then evaluated whether integrating each of these six signatures with our own could improve predictive performance. However, we found that our signature still outperformed the integrated signatures in response prediction in mUC. We further showed that combining our signature and tumor mutation burden (TMB), a widely recognized genomic biomarker of ICI response, significantly improved prediction performance for atezolizumab-treated patients in the IMvigor210 cohort, compared to using TMB alone. Finally, we identified several prognostic markers within our gene signature that were able to stratify overall survival in patients with mUC. Taken together, our method, LogitDA, identified a robust gene signature for mUC that not only outperformed well-established biomarkers in response prediction, but also complemented them to improve the accuracy of immunotherapy outcome prediction.
Materials and methodsData sets
We obtained gene expression profiles and clinical data of patients with mUC and renal cell carcinoma (RCC), respectively, from previously published studies (2, 6). These datasets include clinical response information for 298 patients enrolled in a single-arm phase II clinical trial evaluating atezolizumab as a first-line treatment for primarily mUC patients (IMvigor210) (16), and 77 RCC patients from a randomized phase II trial of atezolizumab vs. atezolizumab+bevacizumab vs. sunitinib in front line RCC (IMmotion150, NCT01984242) (17). To see which of the four classifiers studied can make consistent response prediction to atezolizumab, we used the IMmotion150 and the IMvigor210 datasets as the training datasets in the cross-validation study.
In addition, we obtained data from a phase I clinical trial of atezolizumab (PCD4989g, NCT01375842) (18) through a data request to Genentech (South San Francisco, USA). This dataset includes whole-transcriptome profiles and clinical information from 94 patients with mUC and was used as an independent test dataset. Details regarding the availability of raw data are provided in the Data Availability section. Each dataset contains over 30,000 transcripts with corresponding expression values in raw counts.
Clinical responses in the IMvigor210, PCD4989g, and IMmotion150 datasets were assessed according to the Response Evaluation Criteria in Solid Tumors version 1.1 (RECIST 1.1) (19), and were classified as complete response (CR), partial response (PR), stable disease (SD), or progressive disease (PD). For downstream analyses, patients were grouped into responders (CR/PR) and non-responders (SD/PD), respectively.
Among the 348 patients in the IMvigor210 trial, clinical response data were available for 298 patients, the majority of whom were diagnosed with mUC, although some also had liver, lung, or other cancers. In total, the IMvigor210 (IMmotion150) cohort included 298 (77) patients, with 68 (15) responders and 230 (62) non-responders. The PCD4989g (mUC) test dataset included 94 patients, with 22 responders and 72 non-responders.
RNA-seq data processing and normalization
Whole transcriptomic profiles of IMvigor210 (PCD4989g) were downloaded in raw read (FASTQ) format; the reads were aligned and quantified to gene-level counts by kallisto (20). After applying DEseq2 to identify differentially expressed (DE) genes, we normalized the counts to log2 (TPM + 1) (21) for each sample of training and test sets.
Standardization of training and test data
To reduce the influence of systematic technical differences between training and test data sets, we applied independent standardization (denoted as ST) to IMvigor210 and PCD4989g (mUC), that scaled each gene to mean 0 and variance 1 separately in training and test data, given the proportions of the two outcomes CR/PR and SD/PD are similar between the two datasets.
Feature selection
For feature selection, we have implemented the following procedures. (1) We identified differentially expressed (DE) genes of responders versus non-responders by DEseq2 (22), excluding non-informative genes (HIST and LOC). (2) Next, we calculated the ratio of between-group to within-group sums of squares (BW-ratio) (23) for the DE genes, where the two groups refer to responders and non-responders. and (3) Finally, we applied unsupervised domain adaptation (DA) (10) to sift genes which passed step (2). Details of the feature selection procedures are in Supplementary Materials.
Training the four classifiers by 5-fold CV
After completing the feature selection procedures, we used the IMvigor210 dataset to train the parameters of four classifiers, logistic regression, lasso regression, random forest (RF), and support vector machine with the RBF kernel (SVM(RBF)), via 5-fold cross-validation. For LogitDA, the parameters p of the top-p genes, αDA of DA and the penalty constant λ of logit regression were optimized by grid search with 5-fold CV and 100 repeats to result in the associated predictors. Specifically, p is evaluated from the union of [15(5)200]) and [210(10)500], where the former denotes the set of numbers from 15 with step size 5 to 200 for IMvigor210, and αDA in [0.2, 0.8] (step size = 0.1). Further details of the training procedures for each classifier are provided in the Supplementary Materials.
Statistical analysis and tools
We used the DEseq2 and dgof packages in R software to identify DE genes between responders and non-responders. A log-rank test was applied to reveal genes that can separate patients into favorable or poor OS within IMvigor210. For the remaining analyses, including logistic ridge regression, support vector machine, random forest classification, and visualization (e.g., volcano plots), we used R software.
ResultsWithin-study cross validations demonstrated that LogitDA consistently predicted responses in atezolizumab-treated patients
After applying the feature selection procedures, we performed 5-fold cross-validation (CV) with 100 repeats on the IMvigor210 training dataset using four classifiers, LogitDA, lasso regression, random forest, and support vector machine with the RBF kernel (SVM(RBF)). To determine which classifier provided the most robust predictions across cancer types in atezolizumab-treated patients, we also performed a 5-fold CV study using the IMmotion150 (mRCC) dataset.
In the IMvigor210 dataset, lasso and LogitDA achieved the highest CV AUCs (accuracy) of 0.78 and 0.77 (0.71 and 0.73), respectively, followed by SVM(RBF) and random forest. AUC (area under the receiver operating characteristic curve) was used as the primary metric for classifier performance evaluation. In the IMmotion150 cohort, LogitDA and SVM(RBF) were the top two predictors, each attaining a cross-validated AUC of 0.95.
Detailed CV results of both datasets, including AUC, accuracy, true positives (TPs), true negatives (TNs), and other metrics, are provided in Supplementary Table S1. Based on the results from both datasets, we selected LogitDA and the 49-gene signature it identified for further investigation.
Our signature outperformed the six well-known signatures in ICI-treated patients with mUC
Let IMvigor-PCD(mUC) denote the setting in which a classifier is trained on the IMvigor210 dataset and its prediction performance evaluated on the PCD4989g(mUC) dataset. In this section, using IMvigor-PCD(mUC), we compared the performance of our gene signature with six state-of-the-art signatures, PD-L1 (IHC), IFN-γ (1), tGE8 (2), T exhaust (12), T inflamed (13), and CD8T (14), as introduced in the Introduction. We first conducted 5-fold cross-validation with 100 repeats for each signature using logistic regression combined with our optimization algorithm to derive the corresponding predictors. In the IMvigor210 dataset, our signature achieved the highest cross-validated AUC of 0.77, followed by IFN-γ (0.70), T inflamed (0.69), tGE8 (0.68), and T exhaust (0.62); both CD8T and PD-L1 yielded cross-validated AUCs below 0.60.
Next, we applied the trained predictors to the independent test set, PCD(mUC). Our transcriptomic signature outperformed the six established signatures, achieving the highest prediction AUC of 0.75, followed by T inflamed (0.70), IFN-γ (0.67), tGE8 (0.65), and the remaining signatures (0.59 and lower), as shown in Table 1. In terms of prediction accuracy, our signature also ranked highest at 0.71, followed by IFN-γ (0.67), with the others performing at 0.64 or lower. Since the PCD(mUC) dataset does not include tumor mutation burden (TMB) information, we were unable to evaluate the predictive power of TMB. Detailed metrics, including prediction AUC, accuracy, true positives (TPs), true negatives (TNs), false positives (FPs), and false negatives (FNs), are provided in Table 1. An overview of the training and testing scheme for LogitDA and logistic regression-based predictors, along with their corresponding prediction AUCs in PCD(mUC), is illustrated in Figure 1.
Our signature for mUC outperformed the six known ones in response prediction to Atezolizumab using IMvigor210-PCD(mUC).
Signatures
Parameter
CV result
Prediction result
λa
AUC (SE)
AUC
accuracy
F1-score
TPs
TNs
FPs
FNs
Ours (49)b
0.14
0.77(0.00c)
0.75
0.71
0.47
12
55
17
10
PD-L1d
1.07
0.57(0.001)
0.55
0.40
0.40
15
15
42
3
tGE8 (8)
0.05
0.68(0.01)
0.65
0.61
0.39
12
45
27
10
IFN-γ (18)e
0.05
0.70(0.01)
0.67
0.67
0.46
13
50
22
9
T inflamed (17)e
0.04
0.69(0.02)
0.70
0.64
0.43
13
47
25
9
CD8T (5)
1.04
0.58(0.01)
0.59
0.61
0.39
12
45
27
10
T exhaust (50)
0.15
0.62(0.02)
0.56
0.55
0.28
8
44
28
14
λ is the penalty constant of logistic ridge regression.
LogitDA with the optimized αDA = 0.2 and λ= 0. 15 resulted in the 49-gene model.
SE equals to “0.00” after rounded to the 3rd digit.
After excluding samples with PD-L1 (IHC) missing, the sample size of IMvigor210 and PCD(mUC) were 297 and 75, respectively.
In the IFN-γ and T inflamed signatures, HLA-E was missing in IMvigor-PCD(mUC).
Our signature outperformed the six integrated signatures in response prediction to atezolizumab in mUC
After demonstrating that our signature outperformed the six state-of-the-art signatures individually, we investigated whether combining any of these signatures with ours could further improve predictive performance. As shown in Table 2, our signature achieved the highest prediction AUC (0.75) and accuracy (0.71) in the mUC test set. Specifically, integrating PD-L1 (IHC) and the CD8T signature with ours resulted in prediction AUCs of 0.69 and 0.73, and accuracies of 0.65 and 0.68, respectively. Interestingly, combining PD-L1 with our signature (denoted as PD-L1+ours) led to a reduction in both AUC and accuracy compared to using our signature alone. Nevertheless, the combination CD8T+ours achieved the highest number of true positives (13) among all integrated signatures, while our signature yielded 12 true positives prediction.
Our signature surpassed the integrated signatures in response prediction to Atezolizumab using IMvigor210-PCD(mUC).
Signatures
Parameter
Cross-validation result
Prediction result
p
λa
AUC (SE)
AUC
accuracy
F1- score
TPs
TNs
FPs
FNs
Ours (49)b
49
0.14
0.77(0.00c)
0.75
0.71
0.47
12
55
17
10
PD-L1+oursd
50
0.17
0.78(0.01)
0.69
0.65
0.35
7
42
15
11
tGE8+ourse
56
0.17
0.78(0.01)
0.71
0.66
0.41
11
51
21
11
IFN-γ+oursf
66
0.19
0.77(0.01)
0.71
0.66
0.41
11
51
21
11
T inflamed+oursf
65
0.18
0.77(0.01)
0.72
0.66
0.43
12
50
22
10
CD8T+ours
54
0.18
0.77(0.01)
0.73
0.68
0.46
13
51
21
9
T exhaust+ours
99
0.27
0.75(0.01)
0.71
0.71
0.34
7
60
12
15
λ is the penalty constant of logistic ridge regression.
LogitDA with the optimized αDA = 0.2 and λ= 0. 15 resulted in the 49-gene model.
SE equaled to “0.00” after rounded to the 3rd digit.
After excluding samples with PD-L1 (IHC) missing, the sample size of IMvigor210 and PCD(mUC) were 297 and 75, respectively.
The tGE8 signature having CXCL9 overlapped with our signature.
Of the IFN-γ and Cristescu signaturees, CXCL9 overlapped with our signature, and HLA-E was missing in the IMvigor210-PCD(mUC) datasets.
Differentially-expressed <italic>GABRA3</italic>, <italic>MAST1</italic>, <italic>CXCL</italic>9, <italic>NUF2</italic>, and <italic>LURAP1</italic> were associated with response to Atezolizumab in patients with mUC
From the volcano plot of our 49-gene signature for IMvigor210-PCD(mUC) (Supplementary Figure S1), we identified several genes significantly associated with response to atezolizumab in patients with mUC. We found that GABRA3, MAST1, CXCL9, and NUF2 were the most over-expressed genes, with log2 fold changes of 1.26, 1.24, 1.12, and 0.90, Benjamini-Hochberg adjusted p-values of 6.0×10-4, 8.8×10-7, 0.0025, and 1.4×10-6, respectively). Conversely, LURAP1 was the most underexpressed gene with a log2 fold change of -0.88 and adjusted P =1.5×10-5.
We checked these four over-expressed genes against existing literature, and found the following. Gene GABRA3 has been associated with TMB and shown to promote antitumor immunity in hepatocellular carcinoma based on multi-omics analysis (24). Recent studies have also reported that manipulating GABAergic signaling could limit anti-tumor immunity (25, 26). Since response to immune checkpoint blockade (ICB), such as PD-L1 inhibition, is known to correlate with the extent of tumor immune infiltration, we further investigated immune-related functions of the identified genes. Several reports indicate that CXCL9 is associated with immune cell infiltration (27) and is required for effective antitumor responses following ICB treatment (28).
In addition, Zheng B et al. demonstrated that NUF2 positively correlates with tumor-infiltrating immune cells, including CD8+ T cells and dendritic cells, in clear renal cell carcinoma (29). Notably, our finding that under expression of LURAP1 is associated with better response to atezolizumab is consistent with previous evidence from TCGA bladder cancer data, where hypermethylation of five CpG sites in LURAP1 (resulting in reduced expression) was linked to improved overall survival (30).
Prognostic biomarkers of overall survival identified for mUC
To identify prognostic markers for overall survival (OS), we conducted log-rank tests on the 49 genes included in the LogitDA predictor for the IMvigor210-PCD(mUC) setting. We identified 18 genes as significant prognostic biomarkers, each with a Benjamini–Hochberg adjusted P value (FDR) < 0.01 (log-rank test; Supplementary Table S2). The top five ranked biomarkers were LYRM1, RFC4, CENPL, SPAG5, and CACYBP, all with adjusted P values < 0.0025. Kaplan–Meier survival curves for these five genes are shown in Figures 2A–E.
Kaplan–Meier plots of prognostic biomarkers for overall survival (OS) in the mUC cohort (IMvigor210). Panels (A–E) show survival curves stratified by gene expression levels of (A)LYRM, (B)RFC4, (C)CENPL, (D)SPAG5, and (E)CACYBP, respectively. For each gene, patients were classified into high (≥ median) and low (< median) expression groups within the cohort.
Kaplan-Meier survival curves for five prognostic markers: LYRMI (A), RFC4 (B), CENPL (C), SPAG5 (D), and CACYBP (E). Each plot shows survival probability over time in months, comparing two groups based on gene expression levels (high versus low), with p-values indicating statistical significance.Pathway analysis and functional roles of some signature genes
To elucidate the biological processes captured by our signature, we performed pathway analyses using Reactome (MSigDB) and KEGG. Reactome analysis revealed significant enrichment of MHC class II antigen presentation (P = 0.013) and aberrant mitotic exit regulation (P = 0.029) among the 29 upregulated genes, with marginal enrichment of adaptive chemokine receptor binding (P = 0.08) and immune system pathways (P = 0.09). The 20 downregulated genes were enriched in the cell cycle pathway (P = 0.02). KEGG analysis further identified enrichment in cell cycle, DNA replication, HTLV-1 infection, and DNA repair pathways (adjusted P = 1.9 × 10-16 - 0.005). Details are provided in Supplementary Table S3.
We next examined the functional roles of high-weight and differentially expressed signature genes. Literature evidence supports their involvement in tumor immunity and immunotherapy response. CXCL9, an IFN-γ-inducible chemokine that recruits CXCR3+ effector T and NK cells, is central to anti-tumor immunity (28, 31). In mUC (IMvigor210), elevated CXCL9/IFNG/GBP5 expression correlates with favorable anti-PD-L1 response (32), consistent with the IFN-γ→CXCL9 axis. NUF2 expression positively correlates with tumor-infiltrating CD8+ T cells and dendritic cells in clear renal cell carcinoma, while POLA2 has been identified as a positive biomarker for PD-L1 blockade response in mUC. LURAP1, an adaptor activating the canonical NF-κB pathway, promotes PD-L1 expression and immune evasion; its underexpression in responders aligns with enhanced ICB efficacy (33). CDCA3/5/8, implicated in immune-related pathways (34), were negatively associated with ICI response in our study, consistent with their reported roles in suppressing CD8+ T-cell infiltration. Additionally, SLC6A1 and GABRA3, differentially expressed between responders and non-responders, participate in the GABAergic pathway, whose aberrant activation has been linked to immune suppression in the tumor microenvironment (25, 26).
Ablation study
In this section, we evaluated the contribution of domain adaptation (DA) to the predictive performance of LogitDA. We trained logistic regression models using features selected from the first two feature selection steps, namely excluding the DA filtering step. All model parameters except for αDA were optimized using 5-fold cross-validation on the IMvigor210 dataset. The resulting predictor, denoted as LogitDA–DA (i.e., without DA), included 150 genes. When applied in the IMvigor-PCD(mUC) setting, this model yielded a prediction AUC of 0.63. In comparison, LogitDA (with DA) achieved a 12% improvement in AUC, demonstrating the benefit of incorporating DA.
Discussion
In this study, we developed a feature selection pipeline and a machine learning–based predictor, LogitDA, for identifying a transcriptomic signature of response to PD-L1 inhibition in mUC. We demonstrated that LogitDA robustly predicted patient response to atezolizumab and identified an effective gene signature associated with treatment outcome. Using LogitDA, we further showed that the identified mUC-specific signature outperformed PD-L1 (IHC) as well as five established tumor microenvironment (TME)-associated signatures in terms of both prediction AUC and accuracy.
Furthermore, we integrated our mUC-specific signature with each of the six established immune-related signatures and found that our signature achieved the highest predictive performance, with a prediction AUC of 0.75, outperforming all integrated combinations. These findings suggest that our method is able to identify an effective transcriptomic signature from ICI-treated patients with mUC, and may provide a useful tool for stratifying patients likely to benefit from atezolizumab therapy.
In this study, we developed LogitDA, a feature selection and machine learning–based predictor, to identify a transcriptomic signature of response to PD-L1 inhibition in metastatic urothelial carcinoma (mUC). LogitDA robustly predicted response to atezolizumab and yielded an mUC-specific gene signature that outperformed PD-L1 (IHC) and five established tumor microenvironment (TME)–associated signatures in both AUC and accuracy. When integrated with these immune-related signatures, our signature achieved the highest predictive performance (AUC = 0.75). These results demonstrate that LogitDA effectively identifies clinically relevant transcriptomic predictors and may aid in stratifying mUC patients likely to benefit from atezolizumab.
TMB is a well-established genomic biomarker of response to immune checkpoint inhibitors (ICIs) across several cancer types, e.g., melanoma. Elevated TMB levels are thought to increase neoantigen load, thereby enhancing T cell infiltration and the efficacy of immunotherapy. Unfortunately, the test dataset we assessed did not comprise TMB levels. Thus, we studied the CV result of TMB alone and the combined TMB+ours signature using IMvigor210. The leave-one-out cross-validation AUCs for TMB alone and TMB+ours for mUC were 0.44 and 0.78, respectively. Importantly, we found that the combined TMB+ours signature correctly reclassified 28 non-responders previously misclassified as responders by TMB alone (R2NR), and correctly reclassified 12 responders from previously predicted non-responders by TMB alone (NR2R) for mUC, as summarized in Supplementary Table S4. These results suggest that integrating our signature with TMB can improve the prediction of ICI response in patients with mUC.
Using CIBERSORT (35), we deconvoluted the bulk RNA-seq data from the IMvigor210 cohort to estimate immune cell composition. Among the inferred cell types, plasma cells and M1 macrophages were significantly positively associated with patient response to atezolizumab (P < 0.05, FDR < 14%; Welch’s t-test). To further explore the immunological relevance of our signature for mUC, we compared our genes against the LM22 reference marker gene expression matrix of CIBERSORT. We found that MAST1 and CXCL9 were overlapping genes, and their overexpression was positively correlated with atezolizumab response (IMvigor210).
We acknowledge that the binary classification of responders (Rs) and non-responders (NRs) may limit the clinical granularity of the LogitDA predictor. To evaluate this, we compared the distributions of signature scores among CR/PR, SD, and PD groups using the Mann-Whitney test. The scores of SD were intermediate between those of CR/PR and PD. However, the difference between SD and PD was not significant (P = 0.85), whereas both CR/PR vs. SD and CR/PR vs. PD comparisons were highly significant (P < 2 × 10-11). These results support the appropriateness of binary classification for patient response in the IMvigor210 cohort.
Our methods can be applied to predict ICI treatment response in patients with other cancer types, provided that gene expression data and response outcomes are available for both training and test sets. Furthermore, these datasets are not overly heterogeneous (see also future research directions). Langfelder et al. demonstrated that LogitDA (36), trained on IMmotion150 (a mRCC cohort), identified a 27-gene signature that achieved a prediction AUC (accuracy) of 0.72 (0.83) in PCD4989g (mRCC) (Supplementary Table S5).
To evaluate the clinical generalizability of our signature, we applied it to four independent, unseen datasets GSE176307 (bladder cancer (BLCA), n = 89), GSE111636 (BLCA, n = 11), GSE91061(Melanoma, n = 49), and IMmotion150 (mRCC, n =77), yielding prediction accuracies of 61%, 73%, 67%, and 62%, respectively (Supplementary Table S6). These results may indicate that the signature captures key immune pathways and genes associated with ICI response. However, its predictive performance in external cohorts may depend on the similarity of their 49-gene expression profiles to those of the PCD4989g (mUC) reference set.
Given their translational potential, the 49 signature genes could be developed into a cost-effective diagnostic microarray to evaluate response to ICI therapy in patients with advanced or metastatic UC prior to treatment initiation. To assess the temporal dynamics of patient response and guide adaptive therapy, we analyzed conditional survival curves by condSURV package in R (Q1-Q3 expression levels; Supplementary Figure S2) for the top five prognostic biomarkers. Responders and non-responders began to diverge at 3–4 months, reached maximal separation at 18–21 months, and plateaued thereafter. These results indicate that biomarker discrimination peaks during the mid-term follow-up, providing the greatest clinical utility for monitoring response and informing therapy adjustment within the first 18–21 months.
Several promising gene signatures predicting ICI response have been uncovered by machine learning approaches similar to ours. Shen et al. analyzed glycolysis-related genes (37). They derived an 18-gene signature that achieved a time-dependent ROC AUC of 0.71 in IMvigor210 at 20 months. Notably, several identified genes were significantly associated with response to anti-PD-(L)1 therapy across IMvigor210 and four cohorts. Boll et al. integrated six independent cohorts, encompassing multi-omic data, immune signatures, and others (38). They derived a random forest model, reaching an AUC of 0.76 in a validation set (n = 205) combining IMvigor210 with other cohorts.
We further evaluated three recently reported ICI-response signatures (39–41) using the IMvigor210 (mUC) cohort. The NCOA3/HSP90α/EZH2/CXCL9 axis identified in colorectal cancer by Liu et al. impaired anti-PD-L1 efficacy via CXCL9 suppression; applying this signature (NCOA3, HSP90AA1, EZH2 up, CXCL9 down) for non-responders yielded an accuracy of 0.73. Ke et al. reported PD-1/CD69 co-expression in effector-memory CD8+ T cells as predictive of TLS formation and ICI benefit; classifying mUC patients with both PDCD1 and CD69 above the median as responders achieved 0.57 accuracy. Xu et al. defined a disulfidptosis-related signature (DFRS) associated with poor survival in several cancers but not in bladder cancer (P = 0.177, Fig. 4 (41)); similarly, DFRS failed to stratify IMvigor210 patients (P = 0.092, log-rank) and achieved 0.54 accuracy when high-DFRS cases were predicted as responders.
Recent studies on immunotherapy response prediction in urothelial carcinoma (UC) have explored biomarkers across genetic, proteomic, and transcriptomic levels from tumor, blood, and urine samples. Major research directions include: (1) Combinatorial biomarkers and computational models, particularly ML-based approaches. Yoshida et al. reported that co-expression of LAG-3 and FGL1 was linked to poor PD-(L)1 response and survival, suggesting potential benefit from combined anti-LAG-3/PD-(L)1 therapy (42). (2) DNA damage repair (DDR) alterations, such as mutations in RB1, ATM, BRCA1/2, and ERCC2, which predict ICI benefit or post-chemotherapy responsiveness (43, 44). In a phase II trial, durvalumab plus olaparib failed to improve PFS in unselected mUC patients (n = 154), but higher response rates were observed in those with homologous recombination repair mutations (45). (3) Integration of bulk and single-cell RNA-seq, and increasingly multi-omics data, to identify ICI-associated gene signatures (46–48).
The current optimization method of LogitDA focuses primarily on maximizing cross-validated AUC within the training set. In future work, we plan to apply LogitDA to other cancer types, where the training and test datasets may exhibit greater heterogeneity than those used in this study. Therefore, an important research direction will be to define and quantify the degree of heterogeneity or “distance” between training and test datasets, and to incorporate this measure into stricter criteria for transferring features across studies. Another direction of future investigation is the development of data augmentation techniques to address the class imbalance problem, as the number of responders to ICI treatments is much lower than that of non-responders. Finally, we aim to integrate single-cell RNA-seq with bulk RNA-seq data in mUC; this may elucidate which cell type proportions, and which genes in which cell types are involved in response to ICIs such as atezolizumab.
Data availability statement
The data analyzed in this study is subject to the following licenses/restrictions: Data for the clinical trial IMvigor210 is available in European Genome-Phenome Archive (EGA) under accession number EGAS00001002556. The dataset is also included in the R data package IMvigor210CoreBiologies for the R, accessible at http://research-pub.gene.com/IMvigor210CoreBiologies/. Data for the clinical trials IMmotion150 and PCD4989g are available under restricted access in European Genome-Phenome Archive with the reference number EGAS00001004343. Raw and clinical data of PCD4989g were accessed via Genentech’s internal Strand pipeline (South San Francisco, USA), and downloaded using pEGA3 (https://github.com/EGA-archive/ega-download-client). Requests to access these datasets should be directed to https://ega-archive.org/studies/EGAS00001002556; https://ega-archive.org/studies/EGAS00001004343.
Ethics statement
The IRB-BM committee of Academia Sinica (AS-IRB02-113170) approved this study.
We are grateful the AE and three reviewers for constructive comments that improved our manuscript. We thank Genentech for sharing the raw and clinical data of PCD4989g, Ming-Yueh Huang, Rajat Butola, and Sz-Tsun Hou for discussions, and computational work during the revision.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
The authors declared that Grace S. Shieh was an Associate Editor of Frontiers in Genetics, at the time of submission. This had no impact on the peer review process and the final decision.
Generative AI statement
The author(s) declare that Generative AI was used in the creation of this manuscript. Generative AI was used to tweak or enhance sentences and to correct word/grammar errors.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2025.1607222/full#supplementary-material
ReferencesAyersMLuncefordJNebozhynMMurphyELobodaAKaufmanDR.
IFN-γ–related mRNA profile predicts clinical response to PD-1 blockade. J Clin Invest. (2017) 127:2930–40. doi: 10.1172/JCI91190, PMID: 28650338MariathasanSTurleySJNicklesDCastiglioniAYuenKWangY.
TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature. (2018) 554:544–8. doi: 10.1038/nature25501, PMID: 29443960NassarAHMouwKWJegedeOShinagareABKimJLiuC-J.
A model combining clinical and genomic factors to predict response to PD-1/PD-L1 blockade in advanced urothelial carcinoma. Br J Cancer. (2020) 122:555–63. doi: 10.1038/s41416-019-0686-0, PMID: 31857723ChenZZhouLLiuLHouYXiongMYangY.
Single-cell RNA sequencing highlights the role of inflammatory cancer-associated fibroblasts in bladder urothelial carcinoma. Nat Commun. (2020) 11:5077. doi: 10.1038/s41467-020-18916-5, PMID: 33033240DeleuzeASaoutJDugayFPeyronnetBMathieuRVerhoestG.
Immunotherapy in renal cell carcinoma: the future is now. Int J Mol Sci. (2020) 21:2532. doi: 10.3390/ijms21072532, PMID: 32260578BanchereauRLengNZillOSokolELiuGPavlickD.
Molecular determinants of response to PD-L1 blockade across tumor types. Nat Commun. (2021) 12:3969. doi: 10.1038/s41467-021-24112-w, PMID: 34172722HavelJJChowellDChanTA.
The evolving landscape of biomarkers for checkpoint inhibitor immunotherapy. Nat Rev Cancer. (2019) 19:133–50. doi: 10.1038/s41568-019-0116-x, PMID: 30755690KongJHaDLeeJKimIParkMImS-H.
Network-based machine learning approach to predict immunotherapy response in cancer patients. Nat Commun. (2022) 13:3703. doi: 10.1038/s41467-022-31535-6, PMID: 35764641GaninYLempitskyV eds. (2015).
Unsupervised domain adaptation by backpropagation, in: International conference on machine learning, . Lille:
PMLR.
MourraguiSLoogMVan De WielMAReindersMJWesselsLF.
PRECISE: a domain adaptation approach to transfer predictors of drug response from pre-clinical models to tumors. Bioinformatics. (2019) 35:i510–i9. doi: 10.1093/bioinformatics/btz372, PMID: 31510654YuanSChenY-CTsaiC-HChenH-WShiehGS.
Feature selection translates drug response predictors from cell lines to patients. Front Genet. (2023) 14:1217414. doi: 10.3389/fgene.2023.1217414, PMID: 37519889JiangPGuSPanDFuJSahuAHuX.
Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med. (2018) 24:1550–8. doi: 10.1038/s41591-018-0136-1, PMID: 30127393CristescuRMoggRAyersMAlbrightAMurphyEYearleyJ.
Pan-tumor genomic biomarkers for PD-1 checkpoint blockade–based immunotherapy. Science. (2018) 362:eaar3593. doi: 10.1126/science.aar3593, PMID: 30309915PowlesTKockxMRodriguez-VidaADuranICrabbSJvan der HeijdenMS.
Clinical efficacy and biomarker analysis of neoadjuvant atezolizumab in operable urothelial carcinoma in the ABACUS trial. Nat Med. (2019) 25:1706–14. doi: 10.1038/s41591-019-0628-7, PMID: 31686036LakatosEWilliamsMJSchenckROCrossWCHousehamJZapataL.
Evolutionary dynamics of neoantigens in growing tumors. Nat Genet. (2020) 52:1057–66. doi: 10.1038/s41588-020-0687-1, PMID: 32929288BalarAVGalskyMDRosenbergJEPowlesTPetrylakDPBellmuntJ.
Atezolizumab as first-line treatment in cisplatin-ineligible patients with locally advanced and metastatic urothelial carcinoma: a single-arm, multicentre, phase 2 trial. Lancet. (2017) 389:67–76. doi: 10.1016/s0140-6736(16)32455-2, PMID: 27939400McDermottDFHuseniMAAtkinsMBMotzerRJRiniBIEscudierB.
Clinical activity and molecular correlates of response to atezolizumab alone or in combination with bevacizumab versus sunitinib in renal cell carcinoma. Nat Med. (2018) 24:749–57. doi: 10.1038/s41591-018-0053-3, PMID: 29867230HerbstRSSoriaJ-CKowanetzMFineGDHamidOGordonMS.
Predictive correlates of response to the anti-PD-L1 antibody MPDL3280A in cancer patients. Nature. (2014) 515:563–7. doi: 10.1038/nature14011, PMID: 25428504SchwartzLHSeymourLLitièreSFordRGwytherSMandrekarS.
RECIST 1.1–Standardisation and disease-specific adaptations: Perspectives from the RECIST Working Group. Eur J Cancer. (2016) 62:138–45. doi: 10.1016/j.ejca.2016.03.082, PMID: 27237360BrayNLPimentelHMelstedPPachterL.
Near-optimal probabilistic RNA-seq quantification. Nat Biotechnol. (2016) 34:525–7. doi: 10.1038/nbt.3519, PMID: 27043002WagnerGPKinKLynchVJ.
Measurement of mRNA abundance using RNA-seq data: RPKM measure is inconsistent among samples. Theory Biosci. (2012) 131:281–5. doi: 10.1007/s12064-012-0162-3, PMID: 22872506LoveMAndersSHuberW.
Differential analysis of count data–the DESeq2 package. Genome Biol. (2014) 15:10–1186. doi: 10.1186/s13059-014-0550-8, PMID: 25516281DudoitSFridlyandJSpeedTP.
Comparison of discrimination methods for the classification of tumors using gene expression data. J Am Stat Assoc. (2002) 97:77–87. doi: 10.1198/016214502753479248YinZYanXWangQDengZTangKCaoZ.
Detecting prognosis risk biomarkers for colon cancer through multi-omics-based prognostic analysis and target regulation simulation modeling. Front Genet. (2020) 11:524. doi: 10.3389/fgene.2020.00524, PMID: 32528533ZhangBVogelzangAMiyajimaMSugiuraYWuYChamotoK.
B cell-derived GABA elicits IL-10+ macrophages to limit anti-tumour immunity. Nature. (2021) 599:471–6. doi: 10.1038/s41586-021-04082-1, PMID: 34732892HuangDWangYThompsonJWYinTAlexanderPBQinD.
Cancer-cell-derived GABA promotes β-catenin-mediated tumour growth and immunosuppression. Nat Cell Biol. (2022) 24:230–41. doi: 10.1038/s41556-021-00820-9, PMID: 35145222LiangY-kDengZ-kChenM-tQiuS-qXiaoY-sQiY-z.
CXCL9 is a potential biomarker of immune infiltration associated with favorable prognosis in ER-negative breast cancer. Front Oncol. (2021) 11:710286. doi: 10.3389/fonc.2021.710286, PMID: 34527583HouseIGSavasPLaiJChenAXOliverAJTeoZL.
Macrophage-derived CXCL9 and CXCL10 are required for antitumor immune responses following immune checkpoint blockade. Clin Cancer Res. (2020) 26:487–504. doi: 10.1158/1078-0432.CCR-19-1868, PMID: 31636098ZhengBWangSYuanXZhangJShenZGeC.
NUF2 is correlated with a poor prognosis and immune infiltration in clear cell renal cell carcinoma. BMC Urol. (2023) 23:82. doi: 10.1186/s12894-023-01258-x, PMID: 37138262WangLShiJHuangYLiuSZhangJDingH.
A six-gene prognostic model predicts overall survival in bladder cancer patients. Cancer Cell Int. (2019) 19:1–15. doi: 10.1186/s12935-019-0950-7, PMID: 31516386TokunagaRZhangWNaseemMPucciniABergerMDSoniS.
CXCL9, CXCL10, CXCL11/CXCR3 axis for immune activation–a target for novel cancer therapy. Cancer Treat Rev. (2018) 63:40–7. doi: 10.1016/j.ctrv.2017.11.007, PMID: 29207310WangJZhouYYangYWuZ.
The expression panel of CXCL9, GBP5, and IFNG is a potential pan-cancer biomarker to predict immunotherapy response. Am J Trans Res. (2023) 15:3960., PMID: 37434847BetzlerACTheodorakiM-NSchulerPJDöscherJLabanSHoffmannTK.
NF-κB and its role in checkpoint control. Int J Mol Sci. (2020) 21:3949. doi: 10.3390/ijms21113949, PMID: 32486375JiangDLiYCaoJShengLZhuXXuM.
Cell division cycle-associated genes are potential immune regulators in nasopharyngeal carcinoma. Front Oncol. (2022) 12:779175. doi: 10.3389/fonc.2022.779175, PMID: 35237510NewmanAMLiuCLGreenMRGentlesAJFengWXuY.
Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. (2015) 12:453–7. doi: 10.1038/nmeth.3337, PMID: 25822800LangfelderPLinE-TTsaiY-TChaT-LShiehGS.
Gene signature for response prediction to immunotherapy in mUC and RCC. In: (Tech Rep No 2025-02) Institute of Statistical Science, Academia Sinica, Taipei, Taiwan. (2025). Available online at: https://sitesstatsinicaedutw/library/2025-tech-reports/.
ShenCSuoYGuoJSuWZhangZYangS.
Development and validation of a glycolysis-associated gene signature for predicting the prognosis, immune landscape, and drug sensitivity in bladder cancer. Front Immunol. (2025) 15:1430583. doi: 10.3389/fimmu.2024.1430583, PMID: 39867879BollLMde OcaSVMCamarenaMECasteloRBellmuntJPerera-BelJ.
Predicting immunotherapy response in advanced bladder cancer: a meta-analysis of six independent cohorts. bioRxiv. (2024) 16:589711. doi: 10.1038/s41467-025-56462-0, PMID: 39979258KeHLiPLiZZengXZhangCLuoS.
Immune profiling of the macroenvironment in colorectal cancer unveils systemic dysfunction and plasticity of immune cells. Clin Trans Med. (2025) 15:e70175. doi: 10.1002/ctm2.70175, PMID: 39934971LiuJSuYZhangCDongHYuRYangX.
NCOA3 impairs the efficacy of anti-PD-L1 therapy via HSP90α/EZH2/CXCL9 axis in colon cancer. Int Immunopharmacol. (2025) 155:114579. doi: 10.1016/j.intimp.2025.114579, PMID: 40215778XuSChenZChenXChuHHuangXChenC.
Interplay of disulfidptosis and the tumor microenvironment across cancers: implications for prognosis and therapeutic responses. BMC Cancer. (2025) 25:1113. doi: 10.1186/s12885-025-14246-1, PMID: 40597807YoshidaTNakamotoTAtsumiNOheCSanoTYasukochiY.
Impact of LAG-3/FGL1 pathway on immune evasive contexture and clinical outcomes in advanced urothelial carcinoma. J Immunother Cancer. (2024) 12:e009358. doi: 10.1136/jitc-2024-009358, PMID: 39043605DamrauerJSBeckabirWKlompJZhouMPlimackERGalskyMD.
Collaborative study from the Bladder Cancer Advocacy Network for the genomic analysis of metastatic urothelial cancer. Nat Commun. (2022) 13:6658. doi: 10.1038/s41467-022-33980-9, PMID: 36333289GalskyMDDaneshmandSIzadmehrSGonzalez-KozlovaEChanKGLewisS.
Gemcitabine and cisplatin plus nivolumab as organ-sparing treatment for muscle-invasive bladder cancer: a phase 2 trial. Nat Med. (2023) 29:2825–34. doi: 10.1038/s41591-023-02568-1, PMID: 37783966RosenbergJEParkSHKozlovVDaoTVCastellanoDLiJ-R.
Durvalumab plus olaparib in previously untreated, platinum-ineligible patients with metastatic urothelial carcinoma: a multicenter, randomized, phase II trial (BAYOU). J Clin Oncol. (2023) 41:43–53. doi: 10.1200/JCO.22.00205, PMID: 35737919ChoMChangHKimJH.
Integration of bulk and single-cell RNA-seq reveals prognostic gene signatures in patients with bladder cancer treated with immune checkpoint inhibitors. Cancer Immunol Immunother. (2024) 74:28. doi: 10.1007/s00262-024-03839-7, PMID: 39708127XiaoZLiuXWangYJiangSFengY.
Comprehensive analysis of single-cell and bulk RNA sequencing reveals postoperative progression markers for non-muscle invasive bladder cancer and predicts responses to immunotherapy. Discov Oncol. (2024) 15:649. doi: 10.1007/s12672-024-01548-2, PMID: 39532830LiuXChenCLiJLiLMaM.
Identification of tumor-specific T cell signature predicting cancer immunotherapy response in bladder cancer by multi-omics analysis and experimental verification. Cancer Cell Int. (2024) 24:255. doi: 10.1186/s12935-024-03447-6, PMID: 39033098
Edited by: Zodwa Dlamini, Pan African Cancer Research Institute (PACRI), South Africa
Reviewed by: Zhuming Lu, Jiangmen Central Hospital, China
Xinyu Cao, University of Illinois Chicago, United States