A Novel Computational Framework for Predicting the Survival of Cancer Patients With PD-1/PD-L1 Checkpoint Blockade Therapy.

Background: Immune checkpoint inhibitors (ICIs) induce durable responses, but only a minority of patients achieve clinical benefits. The development of gene expression profiling of tumor transcriptomes has enabled identifying prognostic gene expression signatures and patient selection with targeted therapies.
Methods: Immune exclusion score (IES) was built by elastic net-penalized Cox proportional hazards (PHs) model in the discovery cohort and validated via four independent cohorts. The survival differences between the two groups were compared using Kaplan-Meier analysis. Both GO and KEGG analyses were performed for functional annotation. CIBERSORTx was also performed to estimate the relative proportion of immune-cell types.
Results: A fifteen-genes immune exclusion score (IES) was developed in the discovery cohort of 65 patients treated with anti-PD-(L)1 therapy. The ROC efficiencies of 1- and 3- year prognosis were 0.842 and 0.82, respectively. Patients with low IES showed a longer PFS (p=0.003) and better response rate (ORR: 43.8% vs 18.2%, p=0.03). We found that patients with low IES enriched with high expression of immune eliminated cell genes, such as CD8+ T cells, CD4+ T cells, NK cells and B cells. IES was positively correlated with other immune exclusion signatures. Furthermore, IES was successfully validated in four independent cohorts (Riaz's SKCM, Liu's SKCM, Nathanson's SKCM and Braun's ccRCC, n = 367). IES was also negatively correlated with T cell-inflamed signature and independent of TMB. Conclusions: This novel IES model encompassing immune-related biomarkers might serve as a promising tool for the prognostic prediction of immunotherapy.

Introduction

Tumor cells acquire numerous genomic alterations, deriving “non-self” neoantigens that the immune system can recognize. Although an immune response is noticed in patients with cancer, this response is usually ineffective to tumor elimination (1–4). One of the reasons is the mechanism of immune escape, including profound local immune suppression, induction of dysfunction, tolerance in T-cell signaling, and evasion of immune destruction by the expression of endogenous “immune checkpoints” that generally lead to immune responses after antigen activation (5). These discoveries have increased cancer understanding and developed immunotherapy treatments such as immune checkpoint inhibitors (ICIs) (6). To date, ICIs therapy like anti-PD-1 has been successful for treating many cancers, particularly malignant melanoma, non-small cell lung cancer, and bladder cancer, among others (7–10).

One challenge of the ICIs immunotherapy is the limited proportion of responders, which leads to the urgent need to find predictive biomarkers to identify responders from non-responders (11). Emerging data suggest that patients overexpressing PD-L1 in tumors by IHC have improved clinical outcomes under anti-PD-1 immunotherapy (12). Although PD-L1 IHC seems predictive in lung cancer, it might not be suitable for many other cancers. The microenvironment of tumors has been recognized as a complex system, as the immune response is affected by many different mechanisms besides PD-L1 (13). Besides, IHC-based detection of PD-L1 as a predictive biomarker is confounded by multiple issues, many still unresolved so far, such as variable detection antibodies and cutoff values and the biomarker’s stability and staining of tumor versus immune cells (12–14).

The development of gene expression profiling within tumors has enabled identifying prognostic gene expression signatures and patient selection (15, 16). Recently reported studies had assessed the association of immune-related gene expression in patients with various solid tumors who received immunotherapy. For instance, a genome-wide analysis of melanoma patients treated with recombinant IL2 revealed a signature predictive of clinical response from pretreatment biopsies (14). Moreover, an IFN-inflammatory immune gene expression signature is associated with both enhanced overall response rates (ORRs) and progression-free survival (PFS) in patients with melanoma who received pembrolizumab, which is subsequently being investigated in other malignancies (17). Other examples include an eight-gene signature reflecting pre-existing immunity, the T-effector/IFN-γ signature, explored in a phase II trial of non–small cell lung carcinoma (NSCLC) (18). Although these studies revealed the intrinsic association between pre-existing immunity and the benefit of ICI therapy, the limitations still exist. On the one hand, these studies failed to consider the functional status of pre-existing immunity, which might affect the outcomes of ICI therapy to a large extent. On the other hand, the selection of those signature genes was mainly based on prior knowledge rather than data exploration, which might lead to insufficient application coverage of these signatures.

Previous research revealed two distinct mechanisms of immune escape in tumor (5; Joyce et al., 2015). One is T cell dysfunction, and the other is T cell exclusion. Approaches that measure immune functional signature based on the gene expression profile were developed to explore the correlation with clinical response of immunotherapy, such as Tumor Immune Dysfunction and Exclusion (TIDE) (19), which identified factors which underlie mechanisms of tumor immune evasion.

So far, most published prognostic gene expression signatures have been explored from the perspective of immune activation and elimination (20, 21). However, as another essential character of the tumor microenvironment, immune-exclusive signatures play a suppression role and are rarely researched. Predicting clinical benefit to ICI therapy requires an understanding of how tumors escape the immune system.

In this study, we evaluated the immune-related gene expression profiles in patients with advanced NSCLC, skin cutaneous melanoma (SKCM), and also head and neck squamous cell carcinoma (HNSCC). We are supposed to find an immune signature that can explain the immune-exclusive statement of tumor samples and can predict response to anti-PD-1 checkpoint inhibitor independently of cancer type.

Materials and Methods

Patients and Datasets

This study is a retrospective analysis of patients with immune checkpoint inhibitor therapies. All of the cohorts involved in this study were collected from public datasets, including Prat’s (n=65), Liu’s (n=121), Riaz’s (n=41), Nathanson’s (n=24) and Braun’s (n=181) (22–26). Progression-free survival (PFS) was defined as the beginning of treatment to the date of disease progression (PD). Patients who had not progressed were censored at the date of their last scan. Objective response rate (ORR) was defined as the percentage of patients with complete response (CR) or partial response (PR). Non-objective response (NOR) was defined as the percentage of patients who failed to reach ORR certification. Durable clinical benefit (DCB) was defined as the percentage of patients who achieved CR or PR or stable disease (SD) lasting > 6 months; non-durable clinical benefit (NDB) was defined as PD or SD that lasted ≤ 6 months.

Immune Exclusion Signature

Firstly, the elastic net-penalized Cox proportional hazards (PHs) model was used to select genes with significant power for predictive value in the discovery cohort. Elastic net, a combination of Ridge and most minor absolute shrinkage and selection operator (LASSO) methods, was applied to select prediction features. The regularization parameter, λ, was specified by 10-fold cross-validation, whereas the L1-L2 trade-off parameter, α, was set to 0.5, with equal Ridge and LASSO penalties. The potential prognostic factors determined by the elastic net-penalized CoxPH regression were subjected to multivariate CoxPH regression analysis to adjust the risk scores of each gene chosen for prognostic clinical parameters. Then, variables exhibiting significance in the adjusted analyzes were entered into a backward, stepwise-elimination Cox regression model. The output calculation formula was: RiskScore = gene expression 1×Coef1+gene expression 2×Coef2+…+gene expression n×Coefn:

In this study, the calculation formula of IES was: IES = -2.20*10-5*exp(CCL5)-1.55*10-5*exp(CCR5)+2.67*10-7*exp(CD46)+2.51*10-5*exp(CXCL6)+2.97*10-6*exp(GPI)-4.40*10-4*exp(GZMM)-1.52*10-3*exp(IL13)-1.83*10-3*exp(IL1RAPL2)+2.10*10-6*exp(ITGB1)-1.51*10-4*exp(KLRK1)-7.47*10-5*exp(NFKB2)-1.66*10-4*exp(PDCD1)+1.12*10-4*exp(PLA2G6)-2.56*10-4*exp(TARP)+5.85*10-5*exp(TNFSF4). The software R package cenROC was applied to calculate the ROC of IES for prognostic classification.

Differentially Expressed Genes

The software R package DESeq2 (V.1.30.1) was used to calculate the fold-change of transcripts and to screen for differentially expressed genes (DEGs) (27) in the RNA-seq data. A fold-change larger than two and an adjusted p-value less than 0.05 were set as the cutoff values for screening significant DEGs. Cluster analysis and heatmap generation were performed by the R package and ComplexHeatmap (V.3.12) (25).

KEGG Pathway, GO and GSEA

All differentially expressed genes were subjected to KEGG term analysis and GO biological processes, including calculation of Benjamini-Hochberg corrected p-values through ToppGene(https://toppgene.cchmc.org/) (28). Gene Set Enrichment Analysis (GSEA) was performed using the GSEA software v.3.0 (Broad Institute, Cambridge, USA) (22).

Estimation of Immune-Cell Type Fractions

Cell-type identification by estimating relative subsets of RNA transcripts (CIBERSORT), which is a deconvolution algorithm that can characterize the cell proportion of complex tissues, based on LM22, a normalized gene expression profiles (GEPs) (23, 24, 29). In this study, CIBERSORT (https://cibersort.stanford.edu/) and leucocyte signature matrix 22 (LM22) were used to quantify the proportions of immune-cell types HNSCC samples from the TCGA data. Normalized gene expression data were analyzed by the CIBERSORT algorithm, running 1000 permutations. The CIBERSORT p-value reflects the statistical significance of the results, and a threshold less than 0.05 is recommended. Finally, samples with CIBERSORT p-values less than 0.05 were included in correlation analyzes between genes and immune-cell types.

Immune Gene Signatures

Twenty-three independent gene signatures tracking different cell types (e.g., CD8 T cells, NK cells, and Macrophage) and microenvironment (e.g., cytolytic and dysfunction signatures) were evaluated. Correlation coefficients between IES and different gene signatures were calculated. An unsupervised analysis was performed to cluster correlation coefficients with similar values together. The correlation coefficients of each signature can be found in .

Statistical Analysis

Categorical variables were evaluated with Fisher’s exact tests. Correlation analysis was assessed by Pearson coefficient. Multivariable Cox proportional hazards models were built with gene expressions as covariables. Stepwise regression was used to determine the most informative variables included in multiple (linear) regression models. ROC analysis was done using the cenROC package in R. Significance of overall survival (OS) and progression-free survival (PFS) was determined via Kaplan-Meier analysis with log-rank analysis. The hazard ratio was calculated by the cox function of the survival package in R. All statistical analysis was performed in the R statistical environment version 3.6.1. All tests were two-tailed and a p-value < 0.05 was considered significant.

Results

Immune−Related Gene Expression in the Prognosis of Immunotherapy

To explore the immune-related genes that are related to the prognosis of immunotherapy, we introduced a cohort of 65 patients with advanced NSCLC (n=35), HNSCC (n=5), and SKCM (n=25) from Prat et al. (26). Patients were treated with anti-PD-1 monotherapy, and the expression profile of 730 immune-related genes on this cohort was collected. We compared the relationship between clinical characteristics and immune-related gene expression.

We found that the expression of 41 genes was associated with clinical survival (p < 0.05, respectively) ( ). This 41-gene cluster includes complement-related proteins such as C3, C6 and C8A, C-C/C-X-C chemokine ligands/receptors such as CCL5, CCR5, CXCL6 and CXCR3, Interleukin protein families such as IL2 and IL13, Tumor necrosis factor family/superfamily such as TNF and TNFSF4, immune cell surface markers such as CD8A and CD46. The expression of some genes, such as IFNA17 and IL2, had a higher risk associated with prognosis, while others affected patients’ survival slightly. In order to find the gene combinations and their coefficients that are most suitable for prognosis prediction, we applied “Coxnet”, an algorithm that fits cox model regularized by an elastic net penalty (30). Penalty maximum likelihood estimation was performed with 1000 bootstrap replicates ( ). The optimal weighting coefficient of each gene was determined by the regularisation parameter lambda using the 1–SE standard ( ). Overall, fifteen genes, CCL5, CCR5, CD46, CXCL6, GPI, GZMM, IL13, IL1RAPL2, ITGB1, KLRK1, NFKB2, PDCD1, PLA2G6, TARP and TNFSF4, were selected out to explore patient prognosis ( ). The calculation formula of risk score was defined as:

Figure 1

Construction of the immune−related gene model in the prognosis of immunotherapy. (A) 1000 bootstrap replicates by lasso Cox regression analysis for variable selection. (B) LASSO coefficients of prognosis genes. Each curve represents a prognosis gene. (C) Multivariate Cox proportional-hazards model of 15 immune-related genes based on Prat cohort. Forest plot of 15 immune-related genes and their association with clinical survival, with hazard ratio values shown on the y-axis, and p-values derived from multivariate CoxPH analysis. (D) Hierarchal clustering analysis of 15 immune-related genes in the Prat cohort. '*',P < 0.05. '**', P < 0.01.

RiskScore = -2.20*10-5*exp(CCL5)-1.55*10-5*exp(CCR5)+2.67*10-7*exp(CD46)+2.51*10-5*exp(CXCL6)+2.97*10-6*exp(GPI)-4.40*10-4*exp(GZMM)-1.52*10-3*exp(IL13)-1.83*10-3*exp(IL1RAPL2)+2.10*10-6*exp(ITGB1)-1.51*10-4*exp(KLRK1)-7.47*10-5*exp(NFKB2)-1.66*10-4*exp(PDCD1)+1.12*10-4*exp(PLA2G6)-2.56*10-4*exp(TARP)+5.85*10-5*exp(TNFSF4)

Among 15 genes, expression of six genes, CD46, CXCL6, GPI, ITGB1, PLA2G6 and TNFSF4, increased the risk of distant recurrence, whereas CCL5, CCR5, GZMM, IL13, IL1RAPL2, KLRK1, NFKB2, PDCD1 and TARP expression had a protective effect against prognosis ( ). In particular, each single gene expression had little contribution to the higher or lower risk of prognosis, such as PDCD1 (PD-1, hazard ratio: 0.99), which indicated the complexity of the tumor immune microenvironment ( ).

Prognosis and Clinical Response Prediction With Gene Expression Profile

Fifteen gene expression profiles according to patient prognosis are presented in . All patients were divided into two risk groups according to predict scores based on the regression equation of 15 gene expression profiles. The cutoff of low and high risk was based on the median value of the predicted score (cutoff = -0.133). The result of hierarchical cluster analysis was similar to the predicted-score grouping ( ). Patients in the low-risk group were enriched in expression genes of cancer-suppressing inflammation such as NFKB2 and CCR5. In contrast, patients in the high-risk group were enriched in expression genes of cancer-promoting inflammation such as ITGB1 and CD46.

To understand the relationship between clinical prognosis and 15 risk genes, we analyzed the correlation between progression-free survival (PFS) and individual gene expression ( ). Patients were divided into two groups according to the median expression level. We found that only IL1RAPL2 expression showed a significant correlation in clinical prognosis ( ). The expression of the remaining 14 genes tended to predict prognosis to various extents, though it did not reach significance. When merging the expression of 15 genes, patients in the high-risk group had 2.48-fold higher risk of death compared with patients in the low-risk group (hazard ratio [HR]: 2.48, 95% confidence interval [CI]: 1.34, 4.59, p=0.003, ; ). We applied R package cenROC to analysis the ROC efficiencies and found that the area under the ROC curve (AUC) of 1- and 3- year prognosis were 0.842 and 0.82, respectively ( ). Similar results were observed by cancer types of non-squamous NSCLC and SKCM ( ). Cancer types of squamous NSCLC and HNSCC did not reach significance due to the small candidate size ( ).

Figure 2

Prognosis and clinical response prediction with gene expression profile. (A) Kaplan–Meier survival curves of PFS in high-risk patients versus low-risk patients based on Prat pan-cancer anti-PD-1 monotherapy cohort. (B) Sensitivity and specificity of the prognosis risk score model were assessed by time-dependent ROC analysis. (C) Violionplot of the distribution of risk score value between patients with ORR and NOR. (D) Barplot of object response rate between the high risk group and the low risk group.

Based on the survival analysis results, we evaluated the predicted score and made a prediction model of clinical response in patients. First, we analyzed the predictive performance of every single gene by ROC curve ( ). The objective response’s largest AUC was 0.637 of PLA2G6 gene. None of the gene expression was significantly related to the response of anti-PD-1 therapy, including PD-1 ( , ). However, the predictive risk score value from 15 immune-related gene expression levels was observed to increase significantly. This combination predictive model had an AUC of 0.731, higher than the AUC of PD-L1 expression (AUC=0.625, ). Patients without clinical response had a higher value of predicted score ( ). In particular, only 6 of 33 patients in the high-risk group had an objective response after anti-PD-1 therapy (ORR rate: 18.2%) compared with 43.8% of patients in the low-risk group ( ). Together, the expression profile of 15 filtered risk genes was correlated with the prognosis and clinical response of patients treated with anti-PD-1 therapy.

Pathway and Gene Ontology Analysis Revealed the Difference in Immune Activities Between Two Risk Groups

To identify the inner differences in the tumor microenvironment between two risk groups divided by the predicted score of 15 genes expression, an unsupervised analysis of 730 immune-related genes and risk classification was performed in . We observed that cluster 2, which mainly consisted of patients with a low-risk score, was enriched with a large number of highly expressed immune eliminated cell genes, such as CD8+ T cells (PRF1, CD8A, CD8B, GZMM and FLT3LG), CD4+ T cells (IL26 and IL17A), NK cells (SPN, BCL2 and NCR1) and B cells (BLK and CD19). Then the analysis of differential expression genes was performed between two risk groups, and the expression of 62 genes (of 730) was identified as statistically altered (p <0.05) ( , ; ). The majority of differentially expressed genes displayed decreased expression in the high-risk group. The greatest downregulation of differential expression gene was MS4A1 (2.64 folds), while the expression of ARG1 and S100A7 in the high-risk group upregulated 2.52 folds and 2.66 folds, respectively.

Figure 3

Pathway and gene ontology analysis between two risk groups. (A) Hierarchal clustering analysis of 730 immune-related genes, marked with some markers expressed in CD8+ T cells (PRF1, CD8A, CD8B, GZMM and FLT3LG), CD4+ T cells (IL26 and IL17A), NK cells (SPN, BCL2 and NCR1) and B cells (BLK and CD19). (B) Differential expression analysis between high-risk patients and low-risk patients in the Prat cohort. “UP” indicates that the gene was significantly up-regulated in the high-risk group while “DOWN” indicates the gene was significantly up-regulated in the low-risk group. (C, D) were GO and KEGG enrichment of DEGs, demonstrating that most are related to immune processes. (E) The difference of immune cell infiltration abundances between high- and low-risk patients. (F) and (G) were GO and KEGG enrichment of 15 immune-related genes. (H) Different expression in immune checkpoints (CD274, PDCD1LG2, CTLA4, LAG3, CD28, CD40, CD80, HAVCR2, TIGIT, and TNFRSF9) between high- and low-risk patients. '*',P < 0.05.'**',P < 0.01.'***',P < 0.001.’****’,P < 0.0001. “ns”, no significance.

To better understand the differential expression genes discovered above, we performed pathway enrichment analysis of all differential expression genes between two risk groups by computing their KEGG term and biological process associations. Our analysis generated a total of 19 KEGG terms with a significant p-value (p < 0.05, Benjamini-Hochberg corrected) ( ). Among these KEGG terms, ‘Cytokine-cytokine receptor interaction’ attracted the highest number of differential expression genes, 16 of which seven were discovered among the top twelve differentially expressed transcripts ( ). After gene ontology analysis, we found ten biological process terms that highly correlated with the differential expression genes (p < 0.001, Benjamini-Hochberg corrected) ( ). Most terms were immune-related receptor activity (i.e., T cell, CCR chemokine and CXCR chemokine) and cytokine activity.

Immune cells are essential components of the tumor microenvironment and closely correlate with immunotherapy responses. We used the CIBERSORT software to assess the abundances of 22 different immune-cell types in 65 patients ( ). Similar to the findings of KEGG pathways, naïve B cells, CD4 resting-memory T cells, and activated natural killer (NK) cells account for more enormous proportions of the infiltrating immune cells in the low-risk group than in the high-risk group ( ).

The 15-Gene Expression Profile Is Correlated With Immune-Excluded Microenvironment Characteristics

Path enrichment analysis was also applied to understand the key pathways and biological processes involved in the 15-gene expression profile. Significant pathways and GO terms with two or more enriched genes were selected. The results showed that all of the selected pathways and GO terms were highly relevant to the tumor immune response, such as cytokine signaling, interleukins signaling and lymphocyte activation ( ). These critical steps of immune response explained the expression of 15 risk genes that could predictively evaluate the response to immunotherapy.

To further explore the association between the predictive risk score and various immune-related signatures, Pearson correlation analysis was performed to calculate the pairwise correlations among 24 signatures in 65 patients ( ). The result of hierarchical clustering revealed that our predictive risk score was positively correlated with M2 macrophage signature, cancer-associated fibroblast (CAF) signature and tumor exclusion signature, while negatively correlated with immune-elimination-related signatures such as cytotoxic T lymphocytes (CTL) signature, MHC-II signature and interferon-gamma (IFN-γ) signature. This observation helps us understand that patients with high predicted risk scores received worse clinical outcomes could be explained by the immune exclusion status in their tumor microenvironment to some extends. We also found that the predictive risk score was positively correlated with the expression of classical immune checkpoints, such as PD1, PD-L1, and CTLA4 ( ). Thus, we termed our 15-gene risk score as “Immune Exclusion Score” (IES).

Validation of IES Score in Multiple ICIs Therapy Cohorts

To validate the robustness and eligibility of the IES, we collected three more cohorts that underwent ICIs treatments. All of the patients in the three cohorts were with advanced melanoma, including Liu cohort (n=121), Riaz cohort (n=41) and Nathanson cohort (n=24) (31–33). IES scores were calculated among patients in three cohorts and applied to predict the clinical response of immunotherapy ( – ). ROC analysis showed that the ROC efficiencies of 1- prognosis were 0.547 and 0.622 in Liu and Nathanson cohort, respectively (1- year prognosis AUC was unevaluable in Riaz cohort). While the ROC efficiencies of 3- year prognosis were 0.622, 0.558 and 0.836 in Riaz, Liu and Nathanson cohort, respectively ( ). Patients were divided into IES high and low groups according to the Youden index and the threshold value were -0.0006, -0.0001 and -0.0015 in Riaz, Liu and Nathanson cohort, respectively. We observed a better survival rate and more extended survival advantage in patients within the IES low group in all of the three cohorts ( ). In the Riaz and Nathanson cohorts, the cutoff value was close to the median value of IES among patients. We found that patients with high IES were correlated with worse clinical response in three cohorts, while the objective response rates of the low IES group among three cohorts were 31.8%, 43.9% and 60%, respectively ( ). These results confirmed the predictive performance of clinical outcomes of IES in immunotherapy.

Figure 4

The predictive efficacy of 15 immune-related genes risk score in three validation cohorts. (A–C) Sensitivity and specificity of the risk score model were assessed in each dataset by time-dependent ROC analysis. (D–F) Overall survival analysis between high- and low-risk groups in each cohorts.

Since the cancer type of all three cohorts used for verification above was melanoma, we introduced another immunotherapy cohort of clear-cell renal cell carcinoma(ccRCC) from Braun et al. (34). The ccRCC cohort contains 181 patients treated with anti-PD-1 therapy. Previously studies revealed that traditional biomarkers of immunotherapy such as PD-L1 and TMB didn’t show the ability to distinguish the clinical outcome in ccRCC and the microenvironment of CD8+ T-cell infiltration was related to poor prognosis of immunotherapy. After calculating the IES of ccRCC patients, survival analysis and ROC analysis were applied to evaluate the performance of IES. Intriguingly, our results indicated that patients with high IES, whose tumor microenvironment had the feature of immune exclusion, were correlated with longer PFS (p=0.02, IES cutoff=-0.0976) ( ; ). This result was consistent with the findings in ccRCC, although different from the majority of understanding in the tumor microenvironment. Previous studies proved that PBRM1 mutation was promoting factor of ccRCC ICIs therapy. Thus, we applied IES in the PBRM1-mut subgroup of the ccRCC cohort ( ). The findings in the datasets above showed that patients with high IES scores showed a worse prognosis of OS in the PBRM1-mut subgroup, which indicated that IES could further filter patients with worse clinical outcomes and prognosis ICIs therapy who acquired PBRM1 mutations.

Comparison Between IES and Other Biomarkers of Immunotherapy

Recently, the expression of genes related to cytolytic immune activity was associated with clinical response to ICIs in certain tumors (35, 36). A previous study discovered a T cell–inflamed 18-gene expression profile (GEP) shown to predict response to anti–PD-1 therapy (37). To compare the performance on clinical response to anti-PD-1 therapy between GEP and IES, T cell–inflamed GEP was assessed in all patients from the Prat cohort, Liu cohort, Riaz cohort and Nathanson cohort ( ). We found that higher T cell–inflamed GEP scores were also positively associated with response and prognosis in Prat cohort, Liu cohort, Riaz cohort and Nathanson cohort, showing that the T cell–activated tumor environment also affects response in addition to IES. However, significance was not demonstrated in the subgroups by cancer types in Prat cohort when using GEP as a predictor. We also found that a higher proportion of responders were enriched in the “low-risk” group by IES, compared with GEP.

We next evaluated the correlation between GEP and IES among those patients. Negative correlations were found between GEP and IES, and 57.1% of shared patients were selected as low risk by both GEP and IES ( ). So IES could be a necessary complement to GEP.

Figure 5

Comparison between IES and other biomarkers. (A) Correlation between T cell–inflamed GEP score and IES score in Prat, Liu, and Riaz cohorts. (B) Correlation between somatic mutation counts or TMB and IES score in Liu and Riaz cohorts.

Furthermore, we evaluated the relationship between TMB and IES in Liu and Riaz cohort. IES score was found independent of TMB in the discovery and validation datasets with correlations of -0.11 and 0.011, respectively ( ). This result indicated that the IES score could be applied independently or jointly with TMB in predicting the response of ICIs therapy.

Discussion

The immune checkpoint inhibitors (ICIs) have made remarkable progress in the clinical treatment of tumors in the past decade (38–40). However, immune escape mechanisms and immune resistance to ICIs therapy have not been well-studied. Here, we discovered and developed an immune exclusion-related 15-gene risk score termed “Immune Exclusion Score” (IES) to predict the clinical response and prognosis to anti-PD-1 therapy. Limited clinical outcomes and prognosis to anti-PD-1 treatment occurred in patients with high IES risk scores in the discovery and validation datasets. Besides, IES was also found as a necessary complement to T cell-inflamed GEP signature and independent to TMB, reflective of the relationship of GEP and TMB to IES. These observations suggest that using the IES immune exclusion biomarker may help identify patients who are responsive to anti–PD-1 therapies and explain the mechanism of immune exclusion in the tumor microenvironment.

In this study, we focused on 15 immune-related gene expressions as a prognostic marker of ICIs therapy. Among these genes, CCL5, CCR5 and PDCD1 are essential regulators of T-cell antigen receptor signaling (41–43). Gene GZMM, KLRK1 and NFKB2 activate and improve the cytolytic activity of T and NK cells (44–46). While the high expression of gene CD46, CXCL6 and GPI inhibit the inflammation of the tumor microenvironment (47–49). Our results revealed that CCL5, CCR5, GZMM, IL13, IL1RAPL2, KLRK1, NFKB2, PDCD1 and TARP expression were strong protective effects against distant recurrence of ICIs. The rest genes, CD46, CXCL6, GPI, ITGB1, PLA2G6 and TNFSF4, were revealed for the negative association with tumor suppression.

The differential gene expression and pathway analysis were evaluated between two risk groups divided by the IES scores. Our results showed that patients with low IES score enriched with highly expressed immune eliminated cell genes of CD8+ T cells, CD4+ T cells, NK cells and B cells. Most of the enriched pathways were related to immune-cell membrane receptor activity and cytokine activity. We also observed that naïve B cells, CD4 resting-memory T cells, and activated natural killer (NK) cells infiltrated more enormous proportions in the low IES group than in the high IES group.

Our study results may explain the low response rates and the limited efficacy of ICIs for patients with high IES scores. The IES was positively correlated with immune exclusion signatures such as M2 macrophage signature and CAF signature. In contrast, IES was negatively correlated with immune-elimination-related signatures such as CTL signature, MHC-II signature and IFN-γ signature. In addition, the IES score was positively correlated with the high expression of immune checkpoints like PD1, PD-L1 and CTLA4, which indicated the status of immune suppression among IES-high patients.

The limitations of this study still exist. Firstly, the development and validation of IES were conducted on four retrospective cohorts. Prospective clinical studies are needed to verify the clinical efficacy of IES as a predictive biomarker for immunotherapy. Secondly, the prognostic model is built based only on gene expression data. A model involving more types of data, especially pathological images might be able to improve the prediction accuracy (50, 51). Thirdly, due to the limited size of patients/cohorts, IES was confirmed in a small number of cancer types, such as NSCLC, melanoma, ccRCC and NPC. Data on more cancer types are needed to prove the broad applicability of IES. Besides, all of the cohorts introduced in this study were treated with ICI monotherapy. Combined treatment approaches with ICIs, such as ICI combined with chemotherapy, demonstrated a superior clinical response in recent trials. A more significant implication will illustrate if IES is successfully validated in the combination of ICIs and chemotherapy.

Conclusions

Our data demonstrate that IES can be used to categorize tumors into different subgroups that exhibit distinct patterns of potentially recognizable biology to enhance clinical response. Although the utility of IES, T cell–inflamed GEP and TMB, as well as other emerging agnostic biomarkers, need further validated for response prediction to various immunotherapy approaches, including combination therapies, these findings provide the possibility for further exploring the utility of these biomarkers as guides for clinical precision cancer immunotherapy.

Data Availability Statement

The original contributions presented in the study are included in the article/ . Further inquiries can be directed to the corresponding author.

Author Contributions

XS, HJ, ZG, and ML contributed to research design. Data analysis was carried out by HJ, JW, JX, and XL. JY and XD developed the algorithm and refined the prediction model. ND, HL, and TG collected clinical information that was used in the study and critically revised the manuscript. XS, HJ, and ND wrote the manuscript. All authors provided feedback and approved the final version.

Conflict of Interest

Authors XS, HJ, JW, HL, JX, XL, JY, TG, XD, and ZG were employed by YuceBioTechnology Co., Ltd., Shenzhen, China.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.