Association of transcriptional changes in TEX with ICB response

To delineate the transcriptional changes of TEX during ICB therapy at single-cell resolution, we analyzed single-cell transcriptomic profiles of 52,960 and 43,401 high-quality T cells between pre- and post-ICB TNBC patients from the Zhang cohort [34], and grouped T cells into 12 subsets using the t-SNE (Fig. 1A and Supplementary Fig. 1). After ICB treatment, a significant difference was observed in the relative abundance of exhausted T cells versus other CD8 + T cells (Chi-square test P < 0.001). Specifically, we observed 6546 exhausted T cells (TEX CXCL13 and TEX IFI16) and 2177 CD8 + T cells in pre-ICB patients, and 1733 exhausted T cells (TEX CXCL13, TEX GZMM, and TEX GZMB) and 11,802 CD8 + T cells in post-ICB patients (Fig. 1B). Upon analyzing bulk-tissue transcriptomic profiles of paired ICB-treated melanoma patients from the Riaz cohort, we detected a significant increase in enrichment scores of TEX-related signatures and pathways after ICB treatment (Fig. 1C). However, further investigation of transcriptional changes of TEX with ICB response illustrated a significant increase in enrichment scores of TEX-related signatures and pathways solely in responders, but not in non-responders after ICB treatment (Fig. 1D). Remarkably, there were no significant differences in the enrichment of TEX-related signature and pathways before ICB treatment between responders and non-responders. However, after ICB treatment, responders exhibited higher enrichment scores of TEX-related signature and pathways compared to non-responders (Fig. 1E). These findings substantiate the link between transcriptional changes of TEX and response to ICB immunotherapy.

Fig. 1: Distinct exhaustion profiles of CD8 + T cell during ICB treatment.

A T-distributed stochastic neighbor embedding (t-SNE) visualization of T-cell clusters (52,960 cells for 14 pre-treated patients and 43,401 cells for 12 post-treated patients) with specific markers, showing the annotation and color nodes for T-cell subtypes in the tumor ecosystem. B The proportion of exhausted T cells and other CD8 + T cells between pre-treated and post-treated BC patients (Chi-square test). C Box plots showing computationally estimated activities of TEX signature and pathways in paired ICB-treated melanoma patients (Two-sided paired t test). D Box plots showing changes in the computationally estimated activity of TEX signature and pathways before and after ICB treatment for responders and non-responders, respectively (two-sided paired t test). E Box plots showing the difference in computationally estimated activities of TEX signature and pathways between responders and non-responders before and after ICB treatment (Wilcoxon test).

Deep-learning identification of a TEX-dependent transcriptional signature associated with ICB response

We developed an ensemble deep-learning computational framework, DeepAKR, to identify transcriptional program underlying TEX and ICB response through imbedding the transfer-learning design with supervised pre-training using the RESTIC labeled tumor samples to learn TEX characterizations in melanoma cohort (GHR cohort), followed by parameters fine-tuning on metastatic urothelial carcinoma (mUC) cohort (Maria cohort) to capture the association between ICB response and TEX-related genes (Fig. 2A). Finally, the DeepAKR identified 16 ICB response-associated TEX genes, referred to as ITGs, including SLAMF7, TBX21, IL2RB, IRF1, CCL25, IDO1, GBP4, SLAMF6, CTLA-4, ICOS, SAMD3, ISG20, TIGIT, PDCD1, TOX, and PSMB9 (Fig. 2B). Based on the expression pattern of these genes, unsupervised hierarchical clustering classified 116 melanoma tumors into four clusters (ITG-C1 to ITG-C4) with decreasing expression trend of 16 ITGs from ITG-C1 to ITG-C4 (Fig. 2C). We observed significant differences in T and B lymphocytes composition of the TIME between the four clusters (Kruskal–Wallis test, P < 0.001) (Fig. 2D). Tumors of ITG-C1 and ITG-C2 had high infiltration abundance of lymphocytes, whereas ITG-C3 and ITG-C4 had low infiltration abundance. Furthermore, we scored five TEX-related gene signatures (Cytotoxic, IFN-γ, TNF, IL-2, and CTL) using ssGSEA and found decreasing tendency from ITG-C1 to ITG-C4 (Kruskal–Wallis test, P < 0.001) (Fig. 2E).

Fig. 2: Deep-learning identification of transcriptional program associated with TEX heterogeneity and ICB response.

A Workflow of an ensemble deep-learning computational framework. B Bar plots showing the transcriptional programs with deep learning in the GHR and Maria cohorts. Venn diagram showing the overleaping transcriptional programs between the GHR and Maria cohorts. C Unsupervised hierarchical clustering heatmap of 116 melanoma tumors using expression pattern of 16 ICB response-associated TEX genes (ITGs). D Box plots showing the monotonic association between the computationally estimated abundance of tumor-infiltrating immune cells and ITG subtypes (Kruskal–Wallis test). E Box plots showing the monotonic association between the computationally estimated activity of immune-related biometrics and ITG subtypes (Kruskal–Wallis test). F, G Histograms showing the percentage of each RECIST archetype (CR/PR/SD/PD) among four PD-1 subgroups (F) and ITG subtypes (G) (Chi-squared test). H Pie charts showing the distribution of each RECIST (CR/PR/SD/PD) archetype in four ITG subtypes from the GHR cohort. I Kaplan–Meier curves comparing OS and PFS among four ITG subtypes (log-rank test).

The GHR melanoma solid tumors displayed four major phenotypes using equipartition based on the expression levels of PD-1, and the remaining new subgroups were assigned to these four categories (Fig. 2F). The quantitative analysis of RECIST scores for melanomas revealed that the PD-1 high subgroup had the lowest proportion of progressive disease (PD) state (31.03%), while the PD-1 low cluster exhibited the highest rate of PD state (58.62%). However, there were no significant differences in RECIST scores among the four subtypes of PD-1 expression (Fig. 2F, Chi-squared test P = 0.297). As for the association of the response percentage with ITG clusters, we observed a significant difference in RECIST distribution among four ITG clusters (Fig. 2G, Chi-squared test P = 0.010). The lowest proportion of PD states was found in the ITG-C1 cluster (22.22%), even lower than that in the PD-1 high subtype (31.03%). Meanwhile, the highest proportion of PD states was observed in the ITG-C4 subgroup (73.33%), which was higher than that in the PD-1 low subgroup (58.62%) (Fig. 2H). Moreover, survival analyses revealed that tumors from different ITG clusters in the GHR cohort exhibited significantly different overall survival (OS) and progression-free survival (PFS) (log-rank P < 0.001), with ITG-C1 showing improved considerably survival and ITG-C4 indicating the poorest survival (Fig. 2I).

TEX-dependent machine-learning predictor of response to ICB immunotherapy

Given the observed association of ITGs with ICB response, we developed a TEX-dependent predictor (MLTIP) of response to ICB immunotherapy by utilizing the XGBoost machine-learning method to integrate the 16 ITGs. We trained the MLTIP using the Riaz cohort and then validated the MLTIP in six independent ICB-treated cohorts. We calculated the MLTIP score for each sample and performed a ROC curves analysis using the MLTIP score to assess the predictive power. As shown in Fig. 3A, the MLTIP exhibited superior performance in predicting response to ICB immunotherapy across different cohorts, with AUCs of 0.901, 0.768, 0.762, 0.822, 0.873, 0.750, and 0.573 for the cohorts of Riaz, Gide, HugoW, GHR, Riaz-paired, Nathanson and Maria, respectively. Moreover, the MLTIP scores significantly differed between responders and non-responders, with responders showing higher scores than non-responders. (Fig. 3B). When dichotomizing MLTIP score predictions into either predicted responder-like or predicted non-responder-like groups, the MLTIP showed superior discriminatory power in distinguishing responders from non-responders, with the accuracy of 87.8%, 73.2%, 76.9%, 79.3%, 81%, 88.9%, and 60.7% for the cohorts of Riaz, Gide, HugoW, GHR, Riaz-paired, Nathanson, and Maria, respectively (Fig. 3C). Furthermore, significant differences in OS between predicted responder-like or non-responder-like groups. Tumors in the responder-like group predicted by the MLTIP had significantly better OS compared to tumors classified as non-responders with HR of 0.090, 0.003, 0.868, 0.042, 0.210, 0.010, and 0.267 for the cohorts of Riaz, Gide, HugoW, GHR, Riaz-paired, Nathanson, and Maria, respectively (Fig. 3D).

Fig. 3: Performance of TEX-derived machine-learning predictor of response to ICB immunotherapy.

A ROC curves and corresponding AUC values of MLTIP score. B Box plots showing the distribution of MLTIP scores between responders and non-responders (Wilcoxon test). C Waterfall plot of MLTIP scores and confusion matrices indicating predicted outcomes generated by MLTIP (Chi-square test). D Kaplan–Meier curves comparing OS between responder-like or non-responder-like groups predicted by the MLTIP.

Comparisons to other well-established markers and signatures

We compared the predictive performance of the MLTIP with immune checkpoints (PD-1, PD-L1, and CTLA-4), TMB burden, and six recently proposed signatures in predicting ICB response. The results of ROC analysis demonstrated that the MLTIP consistently achieved superior predictive performance (average AUC = 0.778) compared to other well-established markers and signatures across different cohorts (Fig. 4A, B and Supplementary Fig. S2).

Fig. 4: Performance comparisons to other well-established markers and signatures.

A ROC curves, corresponding AUC values of MLTIP score, and other well-established markers and signatures. B Box plots showing the difference of AUC value between MLTIP scores and other well-established markers or signatures. C Forest plot visualizing HRs of univariate Cox regression analysis of OS in seven cohorts. The red diamond shows the meta-analysis summary of HRs over seven cohorts.

We further compared the prognostic value of the MLTIP with that of six recently proposed signatures in multiple patient cohorts. We first used the univariate Cox regression analysis to evaluate the association of these signatures and the MLTIP with OS for each cohort. Then, we used a meta-analysis to leverage the multiple cohorts for an overall prognostic evaluation of each signature. As shown in Fig. 4C, although the MLTIP and other three signatures (ICB genes, IFN-γ signature, and CTL signature) showed a significant correlation with OS, the MLTIP (HR = 0.093, 95% CI, 0.031–0.280, P < 0.001) demonstrated better predictive performance in OS when compared with ICB genes (HR = 0.832, 95% CI, 0.740–0.935, P = 0.002), IFN-γ signature (HR = 0.329, 95% CI, 0.144-0.751, P = 0.008) and CTL signature (HR = 0.516, 95% CI, 0.320–0.832, P = 0.007).

Association of the MLTIP with immune milieu and prognosis in pan-cancer

We assessed the infiltration of different immune cell subpopulations using the deconvolution methods and enrichment of immune response-related pathways using the ssGSEA across 30 TCGA cancer types, and found that the MLTIP scores highly positively correlated with immune cell infiltration and immune response pathway activation. More specifically, samples with high MLTIP showed increased expression of immune gene signatures (Fig. 5A). This implies that high MLTIP scores could capture active immune tumor microenvironments across various cancer types. In addition, the MLTIP scores vary among patients with the same tumor and between cancer types (Fig. 5B, C), implying general differences in TEX and tumor immunogenicity between different cases of the same tumor and between different tumor types. We further examined the clinical relevance of the MLTIP in pan-cancer, and found that the MLTIP scores decreased significantly for tumors with advanced stages compared with the early stages (Kruskal–Wallis test, P < 0.001) (Fig. 5D). There was a significant difference in the MLTIP scores across gender (Wilcoxon test, P = 0.024), but no difference for different age groups (Kruskal–Wallis test, P = 0.280) (Fig. 5D). The univariate Cox analysis showed that the MLTIP was significantly associated with OS in nine types of cancers (Fig. 5E). However, the prognostic impact of the MLTIP is dissimilar between cancer types. Tumors with high MLTIP phenotype show a significant survival benefit compared with the low MLTIP phenotype in KICH (HR = 0.122, 95% CI, 0.015–0.976, P = 0.018), ESCA (HR = 0.545, 95% CI, 0.315–0.941, P = 0.027), SKCM (HR = 0.484, 95% CI, 0.352–0.665, P < 0.001), BRCA (HR = 0.577, 95% CI, 0.380–0.876, P = 0.009), and BLCA (HR = 0.550, 95% CI, 0.312-0.967, P = 0.035), while the high MLTIP phenotype is mostly associated with reduced OS in LGG (HR = 1.608, 95% CI, 1.137–2.275, P = 0.007), COAD (HR = 2.331, 95% CI, 1.488–3.653, P < 0.001), READ (HR = 3.268, 95% CI, 1.465–7.292, P = 0.002) and UVM (HR = 6.486, 95% CI, 2.603–16.160, P < 0.001) (Fig. 5F).

Fig. 5: Pan-cancer characterization of the MLTIP with immune milieu and prognosis.

A Heatmap of computationally estimated fraction for immune cells (left panel), immune response pathway (right panel) and other immune-related functions (middle panel). MLTIP scores ordered the scheduling of rows. B Alluvial diagram displaying the association between the TCGA cancer type and MLTIP score. C Dot plot showing the distributions of MLTIP scores across 30 cancer types, sorted by the median MLTIP score (horizontal line) for each cancer type. D Box plots displaying the association of MLTIP score and clinicopathologic features (Stage: stage I, II, III, IV; Gender: female, male; Age: young (<=19), adult (20 ~ 39), pre-aging (40–59), old (>60)) with Wilcoxon test and Kruskal–Wallis test. E Forest plot showing the univariate Cox regression analysis of MLTIP score based on OS in TCGA pan-cancer data with two-sided Wald test. The box is presented as HR, and the vertical bar is shown as ± 95% CIs. Green represented the preventive factor, and red represented the risk factor. F Kaplan–Meier curves comparing OS between low-risk and low-risk groups predicted by the MLTIP.