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

Abstract Immune checkpoint therapy (ICB) has conferred significant and durable clinical benefit to some patients with cancer. However, most patients do not respond to ICB, and reliable biomarkers of ICB response are needed to improve patient stratification. Here, we performed a transcriptome-wide meta-analysis across 1,486 tumors from ICB-treated patients and tumors with expected ICB outcomes based on microsatellite status. Using a robust transcriptome deconvolution approach, we inferred cancer- and stroma-specific gene expression differences and identified cell-type specific features of ICB response across cancer types. Consistent with current knowledge, stromal expression of CXCL9, CXCL13, and IFNG were the top determinants of favorable ICB response. In addition, we identified a group of potential immune-suppressive genes, including FCER1A, associated with poor response to ICB. Strikingly, PD-L1 expression in stromal cells, but not cancer cells, is correlated with ICB response across cancer types. Furthermore, the unbiased transcriptome-wide analysis failed to identify cancer-cell intrinsic expression signatures of ICB response conserved across tumor types, suggesting that cancer cells lack tissue-agnostic transcriptomic features of ICB response. Significance: Our results challenge the prevailing dogma that cancer cells present tissue-agnostic molecular markers that modulate immune activity and ICB response, which has implications on the development of improved ICB diagnostics and treatments.


Introduction
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 FDAapproved biomarkers for ICB are PD-L1 expression (3), microsatellite instability (4), and high tumor mutation burden (TMB) of >10 mutations/Mb.(5) While 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.

Clinical Cohorts Treated with ICB
Data from the following studies were used for discovery: i. Kim and colleagues (24) cohort comprises 45 patients with metastatic gastric cancer treated with PD-1 inhibitor pembrolizumab ii.Mariathasan and colleagues (25)  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).

Gene Expression Data Collection and Normalization
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 2 transformed and upper quantile normalized across all cohorts.

Tumor Purity Estimation
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), Pur-Bayes (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 and ICB Response Classification
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.

Expression Deconvolution
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, p 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.

Permutation Statistic to Identify Differentially Expressed Genes
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 Z 2 perm statistics.A P-value is calculated as: That is, the probability of observing Z 2 obs or larger, given the null hypothesis is true.

Identification of Consensus Differentially Expressed Genes Across Cancer Types
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).

Pathway Enrichment Analysis
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 FCERA, GZMH, FCRL, SIRPG, JAKMIP, and WARS as immune-related genes based on literature search.

Analysis of Immune Cell Expression Data from the Human Protein Atlas
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 2 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.ziphttps://www.proteinatlas.org/download/rna_immune_cell.tsv.zip

Analysis of scRNA-seq Data
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.

Ligand-Receptor Analysis
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 ligandreceptor 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.

Predictive Modeling of ICB Response
The IOSelect immune score is defined as the ratio of the mean expression of the top 10 DEGs upregulated in ICB responders (CXCL, CXCL, IFNG, KLRD, GBPP, CRTAM, FASLG, GBP, CDA, PRF) to the mean expression of all eight DEGs downregulated in ICB responders (FCERA, PID, CHH, NPR, PREX, NPR, FAMB, 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 crossvalidation 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 FCERA are the most stably selected features, and built the logistic regression model on all 461 training samples using these three features.

Construction of ICB Response Prediction Model Using Cancer DEGs
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 (RPLL, PGAM, RPSL, TIMM, COPS, COPS, ZNF, TMEMA, and SMAP).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.

Benchmarking ICB Prediction Model Against Existing Biomarkers
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-+TMB model was trained on our training cohort using logistic regression.The cytolytic score is calculated as the geometric mean of GZMA and PRF expression (50).Vigex score is computed as the mean of z-score scaled expression of 12 genes (CXCL, CXCL, CXCL, IFNG, PRF, ILR, GZMA, GZMB, PDCD, CTLA, CD, FOXP; ref. 51).

Immune Cell Enrichment Analysis
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.

Data and Code Availability
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).

Patient Cohorts Used for Transcriptome-wide Analysis of ICB Response
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 metaanalysis 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.

Tumor Transcriptome Deconvolution to Identify Features of ICB Response in Cancer and Stromal Cells
We performed tumor transcriptome deconvolution to estimate the stromaland 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).

Top Stromal Genes Associated with ICB Response are Enriched in Immune Pathways
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).
The top two DEGs with the greatest overexpression in ICB responders were chemokines CXCL and CXCL.CXCL is responsible for T-cell trafficking through binding the CXCR receptor on T cells and has been associated with increased T-cell infiltration in tumors (55).CXCL is involved in both B-cell and T-cell migration via binding to the CXCR receptor.Production of CXCL 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 CXCL and CXCL 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 CXCL and CXCL are indeed the strongest individual transcriptomic predictors of ICB response independent of tumor type.

Stromal Genes Negatively Associated with ICB Response
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, FCERA showed the strongest and most consistent overexpression in ICB nonresponders.FCERA encodes a subunit of the IgE receptor.While FCERA 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).CHH, 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).

scRNA-seq Confirms Differential Expression of the Stromal Biomarkers in Specific Immune Cell Subsets
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 CXCL and CD were most highly expressed in M1-like TAMs, while CXCL and IFNG were enriched in CD8 + T cells (Supplementary Fig. S4).Among DEGs associated with poor ICB response, FCERA was most highly expressed in dendritic cells, mast cells, and M2-like TAMs.PID and CHH were enriched in myeloid cells and fibroblasts, while NPR, NPR,  PREX, 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.

Immune Cell Subtype Analysis Identifies M1 Macrophages as Strong Correlate of ICB Response
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)(69)(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 FCERA, 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).

Cancer Intrinsic Gene Expression Features of ICB Response Tend to be Tumor Type Specific
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 metaanalysis (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.

Lack of Cancer-intrinsic Immune Checkpoint Regulation Across Cancer Types
Although the transcriptome-wide analysis failed to identify conserved cancerspecific 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 (CD) 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.

FIGURE 4
DEGs between responders and nonresponders in the cancer compartment.A, Heat map of genes consistently differentially expressed across discovery cohorts in the cancer cells.Heat maps are colored by the log 2 fold change of gene expression.P-values of individual cohorts are combined using the Fisher method, and the −log 10 (meta P-value) is presented in the bar plot.Genes with q-value <0.01 and P-value <0.1 at least one ICB cohort are shown.B, Bar plot shows the number of significant genes found in the stroma and cancer compartment at different q-value cutoffs.
Pathway analysis of genes with cancer-specific expression differences in gastric cancer (C) and urothelial cancer (E).Bar plots of deconvoluted cancer expression of DEGs in the top enriched pathway in gastric cancer (D) and urothelial cancer (F).Error bars represent estimated SE.
In addition to immune checkpoints, cytokines also play important roles in tumor immunity (74)(75)(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, CXCL, CCL, CXCL, and CXCL), and the corresponding receptors for three of the cytokines (CCR for CCL, CXCR for CXCL and CXCL) 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.

FIGURE 5
Analysis of ligand and receptor expression of immune modulators.Inferred ligand and receptor expression in the cancer and stroma compartments for immune checkpoint-related proteins (A) and cytokines (B).Heat maps are colored by the signed −log10(P-value) of differential expression between ICB responders and nonresponders (permutation test).Tiles with no differential expression (P-value >0.1) were colored gray.In the consensus heat maps, the median signed log 10 (P-value) were plotted for genes that were differentially expressed in three out of four ICB cohorts.The immune checkpoint pairs were annotated by the likely immune cell types as well as the effect of the checkpoint interaction (activating or inhibitory).

A Multivariate Classifier Predicts ICB Response Across Cancer Types
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 immuneactivating 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)(29)(30)(31)(32)(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 FCERA 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.

Discussion
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, CXCL and CXCL, 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, CXCL 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 CXCL expression.In agreement with recent work, we found that CXCL 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, FCERA and CHH 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 FCERA 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 FCERA 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 CCL and CCR 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), triplenegative 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 (FCERA), 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 knowledgeguided 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 cellintrinsic 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.

FIGURE 1
FIGURE 1 Data summary and schematic of workflow for the analysis of determinants of immunotherapy response.

FIGURE 2
FIGURE 2 DEGs and pathways between responders and nonresponders in the stroma.A, Heat map of genes consistently differentially expressed across discovery cohorts in the stroma.Heat maps are colored by the log 2 fold change of gene expression.P-values of individual cohorts are combined using the Fisher method, and the −log 10 (meta P-value) is presented in the bar plot.Genes with q-value <0.01 and P-value <0.1 at least one ICB cohort are shown.B, Bar plots of deconvoluted stromal expression of the top DEGs.Error bars represent estimated SE. C, Pathway enrichment of DEGs in the stroma.Bar plot shows the −log 10 (P-value) of significantly enriched pathways.Similar pathways are grouped by color.D, Network representation of enriched pathways.Large nodes represent enriched pathways.Pathways with high similarity are fused and shown in the same color.Small nodes represent genes associated with each pathway.

FIGURE 3
FIGURE 3 Validation and characterization of stromal DEGs with scRNA-seq and immune cell deconvolution analysis.A, Dotplot showing the percentage of stromal (nonmalignant) cells expressing each gene for each response category.Color represents the average scaled expression of the gene in stromal cells (PD, progressive disease; SD, stable disease; PR, partial response).B, Cell-type specific expression of the top stromal biomarkers in ICB responders (PR) and nonresponders (PD/SD).C, Heat map of the association between immune cell types and ICB response.Abundance of immune cell types of each tumor is estimated using CIBERSORTx.Heat map is colored by the signed −log 10 (P-value) of the difference in CIBERSORTx scores between responders (MSI) and nonresponders (MSS).Cell types with no significant difference in abundance between ICB responders and nonresponders (P-value >0.1) were colored gray.D, Consensus correlations between the abundance of immune cell types and expression of the top stromal biomarkers.Mean Pearson correlation coefficients across the seven discovery cohorts were plotted for each cell type-gene pair.Positive correlations were colored red, negative correlations were colored blue, and tiles with Pearson correlation coefficient <0.1 were colored gray.

FIGURE 6
FIGURE 6 Predictive modeling of response to immune checkpoint inhibition.A, Schematic of data and methods used in model training and testing.B, Top ranked features of ICB response selected by LASSO logistic regression.C, Comparisons of the performance of the IOselect model and the 3-feature model (TMB+IFNG+FCER1A) with existing biomarkers in six independent test cohorts.The differences between the mean AUCs of IOSelect and other models were computed using a paired t test.**, P-value <0.01; *, P-value <0.05; P-value <0.1.D, Kaplan-Meier curves of PFS for patients with high versus low predicted scores (greater than or less than median); P-values from log-rank test shown.

FIGURE 7
FIGURE 7 Schematic summarizing the key components associated with ICB response.Red arrows indicate gene positively correlated with ICB response.Green arrow indicates gene negatively correlated with ICB response.