Transcriptome Deconvolution Reveals Absence of Cancer Cell Expression Signature in Immune Checkpoint Blockade Response

Immune checkpoint blockade (ICB) therapy is an established class of immunotherapy that induces significant tumor shrinkage and long-term disease control in many cancers. However, only a subset of patients benefit from ICB, with response rates ranging from 0% to 3% in pancreatic cancer (1) to about 40% in melanoma and microsatellite instable (MSI) tumors (2). Currently, the FDA-approved biomarkers for ICB are PD-L1 expression (3), microsatellite instability (4), and high tumor mutation burden (TMB) of >10 mutations/Mb. (5) While PD-L1 expression is an obvious biomarker for PD-1/PD-L1 blockade therapies, it is only weakly predictive of treatment response in many solid cancer types (6). TMB is a measure of tumor immunogenicity, where increased mutation load is associated with higher neoantigen burden and antitumoral immune responses. However, the predictive power of TMB varies significantly by cancer type and it is difficult to define a single TMB threshold across cancer types (7, 8). Likewise, MSI tumors are more immunogenic, as defective mismatch repair results in large numbers of insertions and deletions in the cancer genomes. Although MSI status is approved as a tissue-agnostic biomarker of ICB response, it is limited by most solid tumor types not displaying microsatellite instability (9, 10).

Transcriptional signatures of tumor immune phenotypes (11–13) and tumor microenvironment (TME) subtypes (14, 15) have been proposed to predict for ICB response. However, many previous analyses were limited to annotated cancer- or immune-related genes (11, 12, 14, 15) or to single cancer types (13). Furthermore, analysis of large-scale ICB trial data found that existing transcriptional signatures often have limited predictive accuracy in new patient cohorts (6). Two recent meta-analyses investigated the genomic and transcriptomic predictors of ICB response (16, 17). Litchfield and colleagues (16) evaluated previously described biomarkers of ICB response, and found TMB and CXCL9 expression as the most consistent predictors of ICB response.

Bareche and colleagues performed a comparative analysis of existing gene signatures of ICB response, and proposed a new predictive gene signature (17). However, existing large-scale studies have been restricted to bulk tumor transcriptomic data, providing limited insights into the roles of cancer and stroma cells in the TME during ICB therapy.

Activation of effector immune cells such as T cells and natural killer (NK) cells are dependent on a balance of costimulatory and coinhibitory interactions between these effector cells, antigen-presenting immune cells such as dendritic cells and macrophages, as well as cancer cells in the TME. Dysregulation of these ligand–receptor interactions, or immune checkpoints, leads to immunosuppression in the tumor. However, it is challenging to tease apart these complex interactions between cancer and immune cells using bulk tumor transcriptomics. For example, although inhibitory immune checkpoint ligand PD-L1 is a known biomarker of ICB response, it is unclear whether cancer cells escape immune cell killing by expressing PD-L1, or recruiting PD-L1⁺ immune cells, with studies supporting both scenarios (18). Therefore, learning the cell-type specificities of ICB biomarkers can provide insights on the underlying mechanisms of immune evasion and ICB resistance. While single-cell transcriptomics has been applied to explore immune cell types and molecular mechanisms associated with ICB response, single-cell transcriptomics studies are currently limited to small patient cohorts (19, 20). The small sample size, coupled with high intrapatient and interpatient heterogeneity and intrinsic noise in single-cell data, limits the power of single-cell transcriptomics for systematic biomarker discovery.

To study features of ICB response in cancer cells as well as in their neighboring stromal (nonmalignant) cells, we applied a transcriptome deconvolution technique to estimate stroma-specific and cancer cell–specific expression of individual genes from bulk tumor transcriptomes. We used a robust and validated non-negative least squares (NNLS) approach (21), which have previously been applied across 8,000 The Cancer Genome Atlas (TCGA) tumors from 20 solid cancer types to uncover ligand-receptor cross-talk and metabolic states of cancer and stroma cells in the TME (22, 23). Here, we applied this transcriptome deconvolution technique to study the differential expression signatures of ICB responders and nonresponders in cancer and stroma cells. Using a cohort of 1,486 tumors with clinically annotated ICB response and tumors with expected ICB outcomes based on microsatellite status, our analysis revealed a lack of cancer cell–intrinsic gene expression signatures associated with ICB response across tumor types. In contrast, in stromal cells, we identified multiple immunomodulators and checkpoint genes recurrently associated with ICB response across cohorts and tumor types. Many of these genes have previously been linked with ICB response, and we uncovered a subset of novel and potentially immune-suppressive genes negatively correlated with response. Using these conserved stromal-specific gene signatures, we developed a multivariate classifier of ICB response that demonstrated improved accuracy over existing biomarkers in unseen patient cohorts and tumor types.

Data from the following studies were used for discovery:

• Kim and colleagues (24) cohort comprises 45 patients with metastatic gastric cancer treated with PD-1 inhibitor pembrolizumab

• Mariathasan and colleagues (25) cohort comprises 298 patients with metastatic urothelial cancer treated with PD-L1 inhibitor atezolizumab

• Liu and colleagues (26) comprises 118 patients with metastatic melanoma treated with PD-1 inhibitors nivolumab or pembrolizumab

• Gide and colleagues (27) comprises 73 patients with metastatic melanoma treated with anti-PD-1 alone (nivolumab or pembrolizumab), or with anti-CTLA-4 (nivolumab or pembrolizumab with ipilimumab)

Data from the following studies were used for validation:

• Riaz and colleagues (28) comprises 49 patients with advanced melanoma treated with PD-L1 inhibitor nivolumab.

• Pender and colleagues (29) comprises 87 metastatic tumors from 20 solid cancer types treated with a mix of ICB agents.

• Freeman and colleagues (30) comprises 31 patients with advanced melanoma treated a mix of ICB agents.

• Jung and colleagues (31) comprises 27 patients with non–small cell lung cancer treated PD-1 or PD-L1 inhibitor.

• Nathanson and colleagues (32) comprises 19 patients with melanoma treated CTLA-4 inhibitor ipilimumab.

• Ratovomanana and colleagues (33) comprises 24 patients with metastatic MSI colorectal cancer treated nivolumab alone or in combination with ipilimumab.

RNA sequencing (RNA-seq) data are available for all study cohorts, while exome sequencing data are only available for the Kim and colleagues study (24).

Expression data of TCGA colon adenocarcinoma (COAD), rectal adenocarcinoma (READ), stomach adenocarcinoma (STAD), and endometrial tumors were downloaded from the UCSC Xena repository [https://xenabrowser.net/datapages/?cohort=TCGA%20Pan-Cancer%20(PANCAN)&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A44]. Raw transcriptomic data of Kim and colleagues (24) gastric ICB cohort were downloaded from the European Nucleotide Archive, and processed using the bcbio-nextgen pipeline. Briefly, RNA-seq reads were aligned to the human transcriptome using STAR (34), and transcripts per kilobase million (TPM) counts were quantified using SALMON (35). For all other ICB studies, summarized transcriptomic data were obtained from the original publications. Fragments per kilobase million (FPKM) data from Nathanson and colleagues (32) and count data Ratovomanana and colleagues (33) were converted to TPM. For all other ICB cohorts, TPM data were directly downloaded from the original publications. Genes with 0 expression in more than 50% of the samples were removed. All gene expression data were log₂ transformed and upper quantile normalized across all cohorts.

Purity estimates of TCGA tumors were obtained from a previous publication (22), where consensus purity was estimated from mutation allele frequencies, copy-number variations, and RNA profiles of each tumor. As raw genomic sequencing data were unavailable for all ICB cohorts except for Kim and colleagues (36), tumor purity of samples without genomic data were estimated using PUREE (37), a machine learning method for tumor purity estimation from bulk RNA expression data. For the Kim and colleagues gastric ICB cohort (36), where genomic data were also available, we calculated consensus purity estimates based on both genomic and transcriptomic data.

We first computed purity estimates using two RNA-based methods [PUREE (37) and Estimate (38)], and two DNA-based methods [Sequenza (39), PurBayes (40)]. Then, we performed quantile-quantile normalization on the purity estimates of individual methods, and took the mean normalized purity estimate as the consensus purity.

MSI status of TCGA tumors were obtained from two previous pan-cancer studies of microsatellite instability (9, 10). Epstein–Barr virus (EBV) status of TCGA STAD tumors was obtained from a previous TCGA study on gastrointestinal cancers (41). Clinical annotations of ICB response based on the RECIST were obtained from the original studies. Consistent with recent literature (16, 25), we classified complete responders and partial responders as responders and stable disease and progressive disease as nonresponders.

We modeled bulk tumor expression as the sum of cancer and stroma expression, weighted by cancer purity (22). This can be written as:

where e_tumor,i is the bulk tumor expression for the gene in sample i, ē_tumor and ē_stroma are the average cancer- and stroma-specific expression, p_i denotes cancer purity for sample i, and ɛ is the residual. Given bulk expression (e_tumor) and tumor purity (p), NNLS regression can be used to estimate cancer- (ē_tumor) and stroma-specific (ē_tumor) expression for each gene. The SE of the predicted mean expression of compartment h is calculated as:

Where n is the sample size, is the mean tumor purity over all samples, p_h = 1 for cancer-specific expression, and p_h = 0 for stroma-specific expression; MSE is the mean squared error.

To test differences in cancer (or stroma) expression between two response groups (e.g., ē_{tumor, R} – ē_{tumor, NR}), the above approach is applied to each group separately, and a test statistic is calculated as:

To assess significance under the null hypothesis of no difference in cancer (stroma) expression between the two response groups, a null distribution is generated by permutation, where group labels are first randomly shuffled 10,000 times and the above procedure is applied to each permuted dataset, generating a null distribution of statistics. A P-value is calculated as:

That is, the probability of observing or larger, given the null hypothesis is true.

For each gene and each tumor compartment (caner or stroma), a meta P-value is computed by combining the permutation P-values of all cancer cohorts using Fisher method. We used the Bonferroni correction to adjust for multiple testing to maintain the overall α at 0.01. Genes that are differentially expressed in different directions in at least two cohorts were removed. Finally, to eliminate differentially expressed genes (DEG) that are involved in MSI biology but are not directly associated with ICB response, we require the final list of consensus DEGs to be differentially expressed in at least one ICB cohort (nominal P-value <0.1).

To find out whether the DEGs identified are enriched in specific biological processes, pathway enrichment analysis was performed using ClueGo (42) with default parameters. Briefly, two-sided-hypergeometric statistic was used to compute the enrichment of a given list of DEGs in Gene Ontology (GO) biological processes terms, and the resulting P-values were corrected for multiple testing using the Bonferroni step down approach. To aid interpretation, significant biological terms with high overlap as measured by the kappa coefficient were merged into functional groups. To calculate the enrichment of stroma DEGs in immune-related functions, we defined immune-related genes as genes in the “immune system process” term and its child terms in the GO biological processes dataset. In addition, we manually curated stroma DEGs not labeled as immune genes by GO, and further labeled FCER1A, GZMH, FCRL6, SIRPG, JAKMIP1, and WARS as immune-related genes based on literature search.

From the Human Protein Atlas, we downloaded the normalized gene expression of 81 cell types from 31 single-cell RNA-seq (scRNA-seq) datasets, and the normalized gene expression of 18 immune cell types isolated with cell sorting. Then, we plotted the log₂ transformed expression of the candidate stromal DEGs in individual immune cell types.

Download links for Human Protein Atlas data:

• https://www.proteinatlas.org/download/rna_single_cell_type.tsv.zip

• https://www.proteinatlas.org/download/rna_immune_cell.tsv.zip

Gene-level normalized counts from Bi and colleagues (43) were downloaded from the Single Cell Portal. A total of 34,327 cells from eight renal cell carcinoma tumors were analyzed, of which 17,672 cells from five tumors were from ICB-treated patients. Cells with low library size and high mitochondrial contamination were removed. For each annotated cell type, we calculated the average expression and proportion of expressing cells for each DEG, separately for cells from ICB responders and nonresponders. To generate the doplot in Fig. 3A, we scaled the expression counts of each gene using z-score normalization, and computed the average scaled expression and the prevalence of expression of each gene across all stromal (noncancer) cell types.

We curated 30 immune checkpoint ligand-receptor pairs, their cell type specificity and the likely effect of the checkpoint interaction (activating or inhibitory) from literature review. In addition, we curated 25 cytokine/chemokine ligand-receptor pairs involved in immune regulation. To identify ligand-receptor pairs with coordinated expression changes between ICB responders and nonresponders, we surveyed all ligand-receptor genes with nominal P-value >0.1 from the differential gene expression permutation test.

The IOSelect immune score is defined as the ratio of the mean expression of the top 10 DEGs upregulated in ICB responders (CXCL9, CXCL13, IFNG, KLRD1, GBP1P1, CRTAM, FASLG, GBP5, CD8A, PRF1) to the mean expression of all eight DEGs downregulated in ICB responders (FCER1A, PID1, CH25H, NPR3, PREX2, NPR1, FAM134B, SDPR). To build the IOselect model, we used 461 patients from three discovery cohorts [Mariathasan and colleagues (25), Liu and colleagues (26), Kim and colleagues (24)] with both RNA-seq and TMB as training data.

Whole-exome sequencing FASTQ files of Kim and colleagues (24) gastric ICB cohort were downloaded from the European Nucleotide Archive, and processed using the bcbio-nextgen pipeline. Somatic mutations were called by VarScan (44), MuTect (45), VarDict (46), and FreeBayes (47) using default parameters of the bcbio-nextgen pipeline. The final set of mutations was called using an ensemble mutation caller, SMuRF (48), from the outputs of the four mutation callers in the bcbio-nextgen pipeline. TMB was calculated as the number of nonsynonymous mutations per megabase. For the Mariathasan and colleagues (25) and Liu and colleagues (26) cohorts, TMB was obtained from the original publications. All TMB values were standardized using z-score transformation (mean = 0, SD = 1). We fitted a logistic regression model using IOSelect immune score and TMB as input features.

Next, we examined whether we can build a more compact model with comparable performance. Using the same training data, we applied the least absolute shrinkage and selection operator (LASSO) to select top features from the 59 stroma DEGs and TMB. LASSO penalizes the sum of the absolute values of the regression coefficients, shrinking some coefficients to zero, and thereby identifies a simpler and more interpretable model. The regularization parameter λ was chosen by 10-fold cross-validation such that the error of the selected model was within 1 SD from the minimum error. LASSO regression and cross-validation were performed using the “glmnet” package in R. To select for the most robust features, we bootstrapped 100 samples with 90% of the data in each bootstrap, and selected features that are found by at least 75% of the bootstraps for the final model. We found that TMB, IFNG, and FCER1A are the most stably selected features, and built the logistic regression model on all 461 training samples using these three features.

To evaluate whether the cancer-cell specific DEGs identified are predictive of ICB response in independent cohorts, we fitted a logistic regression model using all nine cancer-intrinsic DEGs (RPL22L1, PGAM1, RPS27L, TIMM21, COPS3, COPS4, ZNF717, TMEM176A, and SMAP1). We used patients from all four discovery ICB cohorts as training data and evaluated the performance of the cancer DEG model on six independent test cohorts described in the previous section.

T-cell inflamed gene-expression profile (GEP) scores were calculated from normalized RNA-seq data using the weights provided by Ayers and colleagues 2018 (11). TIDE scores were calculated from normalized RNA-seq data using the TIDE web platform (49). The CXCL-9+TMB model was trained on our training cohort using logistic regression. The cytolytic score is calculated as the geometric mean of GZMA and PRF1 expression (50). Vigex score is computed as the mean of z-score scaled expression of 12 genes (CXCL9, CXCL10, CXCL11, IFNG, PRF1, IL7R, GZMA, GZMB, PDCD1, CTLA4, CD274, FOXP3; ref. 51).

Immune deconvolution was performed using CIBERSORTx (with the LM22 signature matrix; ref. 52) to estimate the absolute abundance of immune cell subsets for each tumor. Immune cell subsets that are differentially infiltrated between responders and nonresponders were identified using the Wilcoxon rank-sum test.

The code used to generate figures in the article is provided in Supplementary Data S1. The transcriptomic and clinical data of the ICB cohorts were obtained from the original publications (Supplementary Table S3).

To systematically explore biomarkers of ICB response transcriptome-wide, we obtained RNA-seq data for 534 pretreatment tumors from ICB-treated patients, comprising four studies of three tumor types [gastric cancer (24), urothelial cancer (25), and melanoma (26, 27); Fig. 1; Supplementary Table S1, and Materials and Methods]. Although microsatellite instability is not a direct measure of ICB efficacy, MSI tumors are highly enriched for ICB responders. Depending on the cancer type, about 30%–50% of MSI tumors respond to ICB, compared with <10% for microsatellite stable (MSS) tumors (2, 53, 54). Adopting a semisupervised learning approach, we considered MSI tumors as ICB responder-enriched, and MSS tumors as nonresponder-enriched, in our meta-analysis for discovery of putative biomarkers of ICB response. We obtained the transcriptomic profiles of 952 tumors from three tumor types in TCGA with the highest frequency of MSI tumors (colorectal cancer, gastric cancer, and endometrial cancer, Fig. 1). As EBV-positive gastric tumors are also enriched for ICB responders (24, 54), we grouped MSI and EBV tumors together as the responder-enriched group for gastric cancer. In total, the discovery cohort comprised of 1,486 tumors from five cancer types.

We performed tumor transcriptome deconvolution to estimate the stromal- and cancer-cell expression of each gene in tumors from responders (MSI) and nonresponders (MSS; Fig. 1). Briefly, we first used a consensus approach to infer tumor purities based on transcriptomic and genomic data (where available, Materials and Methods). For each ICB response group, we estimated gene expression levels in stromal and cancer cells using a NNLS regression approach (22, 23). Finally, we identified cancer/stromal cell DEGs between ICB responders and nonresponders using a permutation-based statistic (Fig. 1; Supplementary Fig. S1). To evaluate the robustness of the deconvolution approach in the ICB cohorts, we confirmed that known lineage-specific genes showed the expected expression differences between cancer and stroma compartments (Supplementary Fig. S2).

From this unbiased transcriptome-wide analysis, we identified 59 genes with differential expression in stromal cells across all cohorts and tumor types (Fig. 2A and B; Supplementary Fig. S3A and S3C, Materials and Methods). To ensure we identify genes that are directly associated with ICB response instead of MSI biology-related genes, we only considered genes that are differentially expressed in at least one ICB cohort. The stroma DEGs identified tended to show consistent differential expression across both MSI/MSS cohorts as well as the four ICB cohorts (Fig. 2A). Among these DEGs, 51 (86%) had immune-related functions as defined by the GO (see Materials and Methods), significantly higher than random expectation (P < 2.2e-16, Fisher exact test). Stromal DEGs include genes involved in immune checkpoint signaling (e.g., CD274, LAG3, CTLA4), IFNγ signaling (e.g., IFNG, IRF1, STAT1), T-cell effector function (e.g., PRF1, GZMA, GZMB), and NK cell–mediated cytotoxicity (KLRD1, KLRC2). Next, we performed pathway analysis to quantify the enrichment of specific biological functions among the stromal DEGs associated with ICB response. We identified 31 enriched GO biological process terms, all of which were immune related (Fig. 2C). These 31 terms could be further grouped into four clusters based on similarity (see Materials and Methods): T-cell activation, chemokine signaling, IFNγ response, and NK cell–mediated cytotoxicity (Fig. 2D).

The top two DEGs with the greatest overexpression in ICB responders were chemokines CXCL9 and CXCL13. CXCL9 is responsible for T-cell trafficking through binding the CXCR3 receptor on T cells and has been associated with increased T-cell infiltration in tumors (55). CXCL13 is involved in both B-cell and T-cell migration via binding to the CXCR5 receptor. Production of CXCL13 by T follicular helper (Tfh) cells is essential for tertiary lymphoid structure (TLS) formation at tumor sites as well as germinal center B-cell activation (56), which has been linked with antitumoral response to ICB (57). Consistent with our meta-analysis, recent studies have also proposed CXCL9 and CXCL13 as key predictors of improved ICB response (58, 59). Remarkably, these two chemokines were also highlighted by a recent meta-analysis that evaluated existing biomarkers of ICB response (16), suggesting that CXCL9 and CXCL13 are indeed the strongest individual transcriptomic predictors of ICB response independent of tumor type.

Although high immune infiltration is generally associated with better ICB response (11, 50, 60), our analysis uncovered a group of immune-related genes consistently associated with poor ICB outcomes across tumor types. Among these genes, FCER1A showed the strongest and most consistent overexpression in ICB nonresponders. FCER1A encodes a subunit of the IgE receptor. While FCER1A is an established marker of basophils, mast cells, and dendritic cells, it is also expressed on immunosuppressive M2 macrophages (61) and tumor-associated macrophages (TAM; ref. 62). CH25H, a gene encoding enzyme cholesterol 25-hydroxylase that catalyzes the formation of 25-hydroxycholesterol (25-HC), was also upregulated in tumors of ICB nonresponders. Interestingly, 25-HC is produced by macrophages in response to type-I IFN signaling and has numerous immunologic effects including B-cell chemotaxis, macrophage differentiation, and regulation of inflammatory response (63). We also observed two natriuretic peptide receptors, NPR1 and NPR3, among genes upregulated in ICB nonresponders. Natriuretic peptides are vasoactive peptides that have been shown to regulate vascular remodeling and promote angiogenesis (64).

Using a scRNA-seq dataset of patients with renal cell cancer treated with ICB (43), we investigated the expression of the stromal DEGs in specific stromal cell subsets and their association with ICB response. First, we confirmed that the stromal DEGs were differentially expressed between ICB responders and nonresponders in stromal (noncancerous) cells (Fig. 3A). Next, we examined the expression distribution of the stromal DEGs among individual cell types. We found that CXCL9 and CD274 were most highly expressed in M1-like TAMs, while CXCL13 and IFNG were enriched in CD8⁺ T cells (Supplementary Fig. S4). Among DEGs associated with poor ICB response, FCER1A was most highly expressed in dendritic cells, mast cells, and M2-like TAMs. PID1 and CH25H were enriched in myeloid cells and fibroblasts, while NPR3, NPR1, PREX2, and SDPR were specifically expressed by endothelial cells (Supplementary Fig. S4). Thus, the stroma DEGs likely encompasses multiple signatures of the TME, including T-cell infiltration, macrophage polarization, angiogenesis as well as fibroblast abundance. Using orthogonal expression data of immune cells subtypes from flow cytometry and scRNA-seq experiments from the Human Protein Atlas, we confirmed the cell-type specific expression of the top stromal DEGs (Supplementary Fig. S5, Materials and Methods). Furthermore, we found that the differential gene expression signal between ICB responders and nonresponders was associated with increased abundance of specific immune cell subsets. For example, the CXCL9+ M1-like TAMs and CXCL13+ T cells were more prevalent in ICB responders than in nonresponders. In contrast, FCER1A+ dendritic cells and FCER1A+ M2-like TAMs were more prevalent in nonresponders (Fig. 3B). Overall, these results suggest that the stromal DEGs represent specific components of the TME that associate with tumor response to ICB.

To further characterize immune cell populations associated with ICB response, we estimated the immune cell composition of each tumor using CIBERSORTx (52) and compared the abundance of immune cell subsets between ICB responders and nonresponders. M1 macrophages, Tfh, and activated CD4 memory T cells had higher abundance in ICB responders as compared with nonresponders (Fig. 3C; Supplementary Fig. S6). M1 macrophages are known to be proinflammatory and antitumoral (65). A recent study found M1 macrophage infiltration to be predictive to ICB response in urothelial cancer (66). Our results further demonstrate that M1 macrophage infiltration could be a general feature of ICB response across multiple cancer types. Tfh cells are specialized CD4⁺ T cells that aid the formation of germinal centers in TLSs (67), interacting with B cells to activate antibody responses (68) and enhancing CD8⁺ T-cell effector functions (69). While the presence of Tfh cells has been associated with favorable survival in multiple cancer types (57, 68–70), their putative role in ICB response has not been studied. Our analysis also demonstrated that while resting CD4 memory T cells were enriched in nonresponders, activated CD4 memory T cells were enriched in ICB responders. This suggests that preexisting CD4⁺ T-cell immunity is important for effective antitumoral response upon ICB treatment.

Next, we examined associations between the top DEGs and infiltration of specific immune cell subsets. We found that the top gene expression markers of good ICB response were generally associated with abundance of M1 macrophages, CD4⁺/CD8⁺ T lymphocytes, and activated NK cells across the discovery cohorts (Fig. 3D; Supplementary Fig. S7; Supplementary Table S2). In contrast, the DEGs associated with poor ICB response, including FCER1A, were correlated with resting state of mast cells, CD4⁺ memory T cells, and dendritic cells, suggesting that these genes are likely markers of a quiescent or immune-suppressive microenvironment (Fig. 3D; Supplementary Fig. S7; Supplementary Table S2).

Previous studies have proposed multiple mechanisms by which cancer cells may directly suppress tumor immunity and impede ICB therapy (18, 71). Surprisingly, our meta-analysis across tumor types revealed a paucity of DEGs associated with ICB response in cancer cells (Fig. 4A and B). In cancer cells, we identified nine genes differentially expressed between ICB responders and ICB nonresponders across tumor types (Fig. 4A; Supplementary Fig. S3B and S3D). However, we observed that the differential expression signal for these nine genes was mostly driven by the MSI ICB responder–enriched cohorts. As compared with stromal DEGs (Fig. 2A), the signal for the cancer DEGs was less consistent across the distinct ICB cohorts (Fig. 4A). Therefore, we hypothesize that the cancer DEGs identified by the meta-analysis are likely related to biology of MSI tumors, rather than general mechanisms of ICB response. Indeed, the nine cancer DEGs were not predictive of ICB response when tested on tumor transcriptomic data from six independent ICB studies (Supplementary Fig. S8). We also confirmed that the paucity of cancer cell DEGs, as compared with stroma cell DEGs, persisted regardless of the significance cutoff in our meta-analysis (Fig. 4B), and we were not able to identify functional enrichment or relationships among cancer DEGs.

Next, we conducted the same DEG analysis within individual tumor types to explore the potential for tissue-dependent processes involved in ICB response. For melanoma, we only identified 22 cancer DEGs shared between the two melanoma cohorts, and these genes were not enriched in specific pathways. In gastric cancer, we identified 190 cancer ICB response DEGs across the gastric ICB cohort and TCGA cohort. Interestingly, genes associated with poor ICB response were enriched for angiogenesis-related functions (Fig. 4C), with proangiogenic genes ENG, F2R, and PDGFRB overexpressed in ICB nonresponders (MSI/EBV) compared with responders (MSS; Fig. 4D). These data are consistent with the understanding that angiogenesis can be associated with immune suppression (72), suggesting that a proangiogenic environment may impede ICB response in gastric cancer. In urothelial cancer, we found that cancer DEGs were enriched for neuron-related functional terms (Fig. 4E), with neuronal genes overexpressed in ICB responders compared with ICB nonresponders (Fig. 4F). This observation is concordant with a previous study linking the neuronal expression subtype (73) of urothelial cancer with poor overall survival but high response rate to ICB (36).

Overall, while these results indicate an overall paucity of tissue-agnostic, cancer cell–intrinsic features of ICB response on the transcriptomic level, they also highlight the potential for tissue and tumor-type specific expression signatures in cancer cells that may predict ICB response.

Although the transcriptome-wide analysis failed to identify conserved cancer-specific DEGs of ICB response, we further explored the expression of immune checkpoint ligands and receptors in cancer and immune cells. We examined the expression patterns of 30 known immune checkpoint ligand-receptor pairs in both the cancer and stromal components of the TME (Materials and Methods), to identify putative ligand-receptor pairs with coordinated expression changes between ICB responders and nonresponders. We found stromal overexpression of 13 checkpoint receptors and two checkpoint ligands to be associated with ICB response across cancer types (Fig. 5A). Checkpoint ligands PD-L1 and PD-L2 and their receptor PD-1 were coordinately upregulated in the stroma of ICB responders compared with ICB nonresponders. PD-L1 (CD274) and PD-1 were significantly upregulated in the stroma of ICB responders in six out of seven cohorts studied (P-value <0.05, permutation test), while PD-L2 was upregulated in five out of the seven cohorts (Supplementary Fig. S9). Interestingly, no checkpoint ligands were consistently differentially expressed and associated with ICB response in cancer cells across tumor types.

In addition to immune checkpoints, cytokines also play important roles in tumor immunity (74–76). We examined the expression of 13 pairs of cytokines and cytokine receptors known to be involved in tumor immunity (74, 76), aiming to identify cytokine signals associated with ICB response (Fig. 5B). Similar to our analysis for immune checkpoints, we found that cancer-intrinsic cytokine expression was not associated with ICB response across tumor types. However, six cytokines (IFNG, FASLG, CXCL13, CCL5, CXCL10, and CXCL9), and the corresponding receptors for three of the cytokines (CCR5 for CCL5, CXCR3 for CXCL9 and CXCL10) were significantly upregulated in the stroma of ICB responders, consistent with previous observations (11, 26). Overall, our analysis finds no evidence for a universal mechanism by which cancer cells inhibit immune cells via checkpoints interactions or chemokine/cytokine signaling.

Next, we explored whether we could build an accurate predictor of ICB response based on the consensus DEGs identified earlier (Fig. 6A). We focused on the 59 stroma-specific DEGs (Fig. 2A) as they were more consistent across tumor types, in comparison with cancer DEGs, and therefore more likely to be generalize to unseen cancer types. To capture the balance of immune-activating and suppressive signals in the TME, we developed the IOSelect immune score, capturing the relative expression difference of response and resistance-associated DEGs (Materials and Methods). While a stromal transcriptomic signature reflects the state and favorability of the TME, TMB is a proxy of the tumor immunogenicity. Because the stromal DEGs and TMB likely represent different facets of tumor immunity, we hypothesized that these could be orthogonal predictors of ICB response. Therefore, we built a multivariate classifier, IOSelect, using both the IOSelect immune score and TMB as input features. The IOSelect classifier was trained on 461 samples from the ICB discovery cohorts with both transcriptomic and genomic data. We evaluated the model performance on six independent test cohorts treated with ICB (28–33), comprising 232 patients with melanoma, colorectal, lung, and other cancers (Materials and Methods). We calculated the area under the receiver operating curve (AUC) of the model, and benchmarked it against established biomarkers of ICB response, namely CXCL9+TMB (16), T-cell GEP (11), VIGex (51), cytolytic score (50), TIDE (13, 49), TMB alone (2, 8), and PD-L1 expression (3). The IOSelect classifier achieved a mean AUC of 0.72, significantly outperforming existing biomarkers with mean AUCs ranging from 0.60 (TMB alone) to 0.68 (T-cell GEP; Fig. 6B).

Next, we evaluated whether a simpler model comprising just the top few predictive features could be sufficient to robustly predict ICB response. We used LASSO logistic regression with 10-fold cross-validation to select the most robust predictors from the set of 59 stromal DEGs and TMB (Fig. 6A). We found TMB, IFNG expression, and FCER1A expression to be the most frequently selected features from 100 bootstrap trials (Fig. 6C). The compact 3-feature model attained a mean AUC of 0.70 across the six test cohorts, which is only marginally lower than the AUC attained by the full IOSelect model (Fig. 6B, P-value = 0.40 by paired t test). Furthermore, the 3-feature classifier is predictive of progression-free survival (PFS) in cohorts with survival data available (Fig. 6D) and to a lesser extent, overall survival (Supplementary Fig. S10). High predicted score is significantly associated with longer PFS in the Gide (27), Liu (26), and Pender (29) cohorts, and is positively associated with PFS in the Ratovomanana (33) and Freeman (30) cohorts although statistical significance is not reached, likely due to the smaller sample sizes of these two cohorts (Fig. 6D). Overall, these results demonstrate that the novel negative biomarkers of ICB response identified from our transcriptome-wide analysis contributes additional predictive power over existing biomarkers.

In this study, we performed a transcriptome-wide meta-analysis of ICB response in 1,486 tumors. These tumor samples included 534 tumors with annotated clinical response to ICB treatment, and 952 tumors from cancer types enriched for MSI and ICB-responding tumors. Using in silico tumor transcriptome deconvolution, we analyzed gene expression signatures of cancer and their neighboring stromal (nonmalignant) cells and identified DEGs associated with ICB response across tumor types. Strikingly, this meta-analysis revealed a paucity of ICB-response DEGs in cancer cells as compared with stromal cells. The few cancer DEGs that we identified by the meta-analysis were more likely related to the biology of MSI tumors and their DNA mismatch repair deficiency rather than being general determinants of ICB response. In contrast, our analysis revealed 59 stroma-specific ICB-response DEGs, which were strongly enriched in immune-related pathways and processes (Fig. 7).

The two chemokine ligands, CXCL9 and CXCL13, constituted the top two stromal genes overexpressed in ICB responders. Underlining the validity of our approach, these two genes have been highlighted for their roles in antitumor immunity and association with favorable ICB response in multiple recent studies (16, 58, 59). Previous studies suggest CXCL13 may play dual roles in ICB response. First, CXCL13 is produced by Tfh cells and involved in B-cell organization in TLS (56), and the presence of TLS can be a positive predictor of ICB response (57, 69). Second, CXCL13 has been reported to be highly expressed by neoantigen-reactive CD8⁺ T cells directly involved in cancer cell killing (77, 78). In line with these findings, we found enrichment of both Tfh and CD8⁺ T cells in tumors of ICB responders compared with nonresponders across cancer types (Figs. 3 and 7). In addition, the abundance of these two cell types correlates with CXCL13 expression. In agreement with recent work, we found that CXCL9 expression is strongly correlated with M1 macrophage infiltration (58, 79), and M1 macrophages were one of the most differentially enriched immune cell types in ICB responders (Fig. 3A). Although myeloid infiltration has traditionally been associated with immune suppression (65), our data are consistent with the hypothesis that M1 macrophage infiltration, and macrophage polarization, could be an important determinant of tumor sensitivity to ICB.

Intriguingly, we also identified a set of immune-related genes that were negatively associated with ICB response. Among these genes, FCER1A and CH25H have previously been reported to be overexpressed in anti-inflammatory M2 macrophages (61, 80), and FCER1A+ TAMs have been reported to promote tumor progression by engaging in a positive feedback loop with tumor-initiating cells (62).

Interestingly, we observed an enrichment of the resting states of CD4 memory T cells, dendritic cells, and mast cells in ICB nonresponders. In contrast, the activated states of these cell types are enriched in the tumors of ICB responders (Fig. 3B). Furthermore, while FCER1A is expressed in antigen-presenting cells like dendritic cells and macrophages, we noted that its expression is 25 times higher in resting dendritic cells compared to activated dendritic cells, and 14 times higher in M2 macrophages compared with M1 macrophages (52). These results suggest that a preexisting immunosuppressive microenvironment, associated with resting immune cells and high FCER1A expression, may prevent effective antitumor immune responses in ICB nonresponders.

Analyses of immune checkpoint and cytokine ligand–receptor interactions revealed coordinated upregulation of several ligand-receptor pairs in the stroma of tumors from ICB responders, including genes encoding PD-L1/PD-L2 and their receptor PD-1, CXCL9/CXCL10 and their receptor CXCR3, and CCL5 and its receptor CCR5 (Figs. 5 and 7). In contrast, we did not find consistent expression changes of checkpoint proteins and cytokines in cancer cells of ICB responders. Among chemokines upregulated in ICB responders, CXCL9, CXCL10, and CXCL13 promote lymphocyte recruitment, while CCL5 is a chemoattractant of both T cells and myeloid cells. CCL5 also directs the recruitment and differentiation of CCR5+ blood monocytes and neutrophils into immunosuppressive TAMs and neutrophils (74, 75). In addition, CCL5 has direct protumoral effects on cancer cells, promoting cancer proliferation and metastasis (74, 75). Preclinical and clinical studies suggested that targeting the CCL5-CCR5 axis could enhance antitumoral response (81, 82). Furthermore, CCR5 was found to be highly expressed by neoantigen-reactive T cells, and overexpressed in ICB responders (16). Here we found both CCL5 and CCR5 to be overexpressed in ICB responders, suggesting that targeting CCL5-CCR5 signaling in conjunction with ICB might further potentiate antitumoral immune response. Indeed, a recent trial found modest clinical benefit of anti-CCR5 and anti-PD-1 combination on patient with heavily pretreated microsatellite stable colorectal cancer, and more trials are underway to evaluate the safety and efficacy of anti-CCR5 and checkpoint inhibitor combinations (83). Currently, the majority of companion diagnostics for PD-1/PD-L1 inhibitors measure PD-L1 expression level on either cancer cells alone, or on both cancer cells and infiltrating immune cells, but seldomly on immune cells only (18, 84). Recent clinical trials found PD-L1 expression on immune cells to be prognostic in patients with head and neck cancer (85), non–small cell lung cancer (86), triple-negative breast cancer (87), and renal cell carcinoma (88) treated with ICB. These data are consistent with our results, suggesting that PD-L1 expression on immune cells could be a more consistent marker of ICB response across tumor types. Using the stromal consensus DEGs identified from the meta-analysis, we developed the IOSelect score, which robustly predicted ICB response in six independent test cohorts, outperforming existing gene expression signatures for ICB response. Furthermore, we proposed a compact 3-feature model that can attain comparable performance with a single positive marker of immune infiltration (IFNG), a single negative marker of immune suppression (FCER1A), in combination with TMB. While ideal pan-cancer biomarkers of ICB response should be robust across patient cohorts, we noticed that the Liu and colleagues melanoma cohort (26) did not share many DEGs with the other cohorts, including other melanoma cohorts. This result highlights significant heterogeneity between ICB cohorts, even among cohorts of the same tumor type. It is unclear to what extent this heterogeneity arise from patient selection, differences in ICB regime, differences in prior treatments, different sample collection protocols, or differences in gene expression assays. This highlights the need for large, consistently selected patient cohorts, with standardized gene expression profiling to facilitate systematic biomarker discovery. Furthermore, combining unbiased transcriptome-wide biomarker discovery with prior knowledge-guided analysis of gene signatures may further improve patient selection for ICB.

As cancer cells are under negative selective pressure from the immune system, several mechanisms of immune escape by cancer cells have been proposed: (i) cancer cells downregulate antigen presentation by deleting HLA alleles (71), (ii) cancer cells deplete neoantigens to avoid detection (89), (iii) oncogenic signaling may upregulate PD-L1 expression and suppress T-cell recruitment (90). However, the prevalence of these immune evasion mechanisms are still unclear. A recent meta-analysis found no association between LOH of HLA alleles and ICB response across ICB cohorts (16). Furthermore, a pan-cancer analysis failed to detect neoantigen depletion signals in most tumor types (91). Similarly, we found a lack of conserved cancer-intrinsic mechanisms of immune evasion on the transcriptomic level, which challenge the prevailing dogma that cancer cell–intrinsic expression of checkpoint ligands can modulate immune activity and predict ICB response across cancer types. These insights have profound implications on biomarker assays and in vitro model systems for development of improved ICB drugs. These results suggest that cancer cells adopt diverse mechanisms of immune evasion, with no universal strategy adopted by cancer cells to avoid immune-mediated killing. The growing number of molecularly characterized ICB cohorts will enable further exploration of the importance and prevalence of the different cancer-intrinsic immune evasion strategies in determining ICB response.

Y.A. Guo reports a patent to Data-Driven Immune Checkpoint Blockade Therapy Response Prediction pending. A.C. Tan reports personal fees from Amgen, Bayer, and Pfizer outside the submitted work. D.S.W. Tan reports personal fees from Takeda, Roche, Merck, DKSH, Bayer, Genmab, and Zymeworks; grants and personal fees from Pfizer, AstraZeneca, and Amgen outside the submitted work. I.B. Tan reports personal fees from Roche, MSD, BMS, and Novartis outside the submitted work. A.J. Skanderup reports a patent to Data-Driven Immune Checkpoint Blockade Therapy Response Prediction pending. No disclosures were reported by the other authors.

Y.A. Guo: Data curation, formal analysis, supervision, funding acquisition, investigation, visualization, methodology, writing-original draft, writing-review and editing. T. Kulshrestha: Investigation, visualization, methodology. M.M. Chang: Investigation, visualization. I. Kassam: Methodology. E. Revkov: Methodology. S. Rizzetto: Investigation. A.C. Tan: Funding acquisition, investigation, writing-review and editing. D.S.W. Tan: Investigation, writing-review and editing. I.B. Tan: Investigation, writing-review and editing. A.J. Skanderup: Conceptualization, funding acquisition, investigation, methodology, writing-original draft, project administration.

This study was supported by the Singapore Ministry of Health's National Medical Research Council under its OF-IRG program (OFIRG18may-0075) and Singapore Agency for Science, Technology and Research (A*STAR) under its Biomedical Engineering Programme (C221118001 and C221318004)