Kun Zheng
Lianchong Gao2†
Xin Zou- 1Department of Urology, Shanghai Jiao Tong University Affiliated Sixth People’s Hospital, Shanghai, China
- 2Yantai Institute, China Agricultural University, Yantai, China
- 3Institute of Clinical Science, Zhongshan Hospital, Fudan University, Shanghai, China
- 4Center for Tumor Diagnosis & Therapy, Jinshan Hospital, Fudan University, Shanghai, China
- 5Department of Pathology, Jinshan Hospital, Fudan University, Shanghai, China
Most patients with clear cell renal cell carcinoma (ccRCC) have an impaired response to immune checkpoint blockade (ICB) therapy. Few biomarkers can predict responsiveness, and there is insufficient evidence to extend them to ccRCC clinical use. To explore subtypes and signatures of immunocytes with good predictive performance for ICB outcomes in the ccRCC context, we reanalyzed two ccRCC single-cell RNA sequencing (scRNA-seq) datasets from patients receiving ICB treatment. A subtype of proliferative CD4+ T cells and regulatory T cells and a subtype of antigen-presenting monocytes that have good predictive capability and are correlated with ICB outcomes were identified. These findings were corroborated in independent ccRCC ICB pretreatment bulk RNA-seq datasets. By incorporating the cluster-specific marker genes of these three immunocyte subtypes, we developed a prediction model, which reached an AUC of 93% for the CheckMate cohort (172 samples). Our study shows that the ICB response prediction model can serve as a valuable clinical decision-making tool for guiding ICB treatment of ccRCC patients.
Introduction
Clear cell renal cell carcinoma (ccRCC) is the most common histological subtype and accounts for more than 70% of all renal cell carcinoma (RCC) (1). The global incidence has shown an upward trend in recent years (2). Large-scale tumor transcriptome analysis revealed that ccRCC was highly infiltrated by immune cells (3), and most notably, its T-cell infiltration was the highest among tumor types within The Cancer Genome Atlas (TCGA) (3). Based on this feature, immune checkpoint blockade (ICB) therapy has significantly improved advanced ccRCC treatment (4, 5). However, a large proportion of ccRCC patients do not respond to this therapy for unclear reasons (6). Such diversity in response to ICB highlights the necessity of identifying predictive biomarkers.
Previous whole-exome sequencing (WES) and transcriptome sequencing of tumors identified several factors associated with ICB outcomes, such as high T/low myeloid cell infiltration (7), high B-cell (8)/CD4 memory T cell abundance (9), elevated PD-L1 expression (10), high tumor mutational burden (TMB) (11), and high similarity/diversity of TCR clonality (12, 13), which are associated with ICB responses. Exhausted/dysfunctional CD8+ T cells (13, 14), anti-inflammatory/M2-like tumor-associated macrophages (TAMs) (13, 14), defects in IFNγ signaling (15) or antigen processing and presentation (16) all contribute to resistance. However, ccRCC is unique compared to other tumors in response to immunotherapy. For example, CD8+ T-cell infiltration (17), tumor mutation burden (11), frameshift insertions and deletions (fsINDELs) (18), and HLA heterozygosity (19) have been described as pro-response effects in other tumor types, whereas none of them have been shown to be associated with ICB response in ccRCC (4, 7, 20). To date, the factors driving ICB resistance remain largely unknown. Effective biomarkers for predicting the ccRCC ICB response in clinical practice are still lacking.
In this study, we reanalyzed two independent ccRCC single-cell RNA sequencing (scRNA-seq) datasets (12, 14) to construct a prediction model. We found that proliferative CD4+ T cells (MKI67+ CD4Ts) and proliferative regulatory T cells (MKI67+ Tregs) can promote ICB resistance, and a subset of antigen-presenting monocytes is related to the ICB response. Based on these findings, we developed an effective ICB response signature, and this signature was validated in several ICB pretreatment bulk RNA-seq datasets. Overall, this signature reached an AUC of 93% for a 172-sample cohort and is more effective than previously published ICB response signatures. These findings extend our understanding of factors associated with the ICB response and provide a potentially powerful response prediction model for ICB clinical treatment.
Results
CD4+ T, Treg, and monocyte cells are associated with the ICB response/resistance in ccRCC
We reanalyzed two independent publicly available ccRCC scRNA-seq datasets [Bi’s dataset (14) and Au’s dataset (12)], which were both from ccRCC patients treated with PD-1/PD-L1 monoclonal antibodies. Then, the observations from single-cell datasets were validated across multiple ccRCC bulk RNA-seq datasets.
Details of all scRNA-seq and bulk RNA-seq datasets used in this study are described in the Methods, and the metadata for all samples are available in Supplementary Table 1. Patients were divided into responders (R), including complete response (CR) and partial response (PR), and nonresponders (NR), including stable disease (SD) and progressive disease (PD), according to the response evaluation criteria in solid tumors (RECIST) (21). Not evaluable (NE) samples from all cohorts were omitted from the RECIST response analysis but remained from the survival analysis.
The overall design of this study is summarized in Figure 1. We used the cell types of Bi’s dataset (n.patients = 7, R = 2, NR = 2, NoICB = 3) and Au’s dataset (n.patients = 2, R = 1, NR = 1) defined by the original studies (Figure S1A). There were ten cell types from Bi’s dataset, i.e., B cells, CD4+ T cells (CD4Ts), CD8+ T cells (CD8Ts), dendritic cells (DCs), monocytes, tumor-associated macrophages (TAMs), natural killer cells (NKs), natural killer T cells (NKTs), regulatory T cells (Tregs) and cancer cells, as well as three cell types from Au’s dataset: CD4Ts, CD8Ts and Tregs. Figure S1B shows the visualization results of the 12 immune cell types (excluding cancer cells) in two-dimensional space via uniform manifold approximation and projection (UMAP).
Figure 1 The overall design of this study. Au and Bi are scRNA-seq datasets, and CheckMate, CM025, CM010, CM009, Javelin 101, GSE67501 and PMC6035749 are bulk RNA-seq datasets. GSEA, gene set enrichment analysis; DE, differentially expressed; ICB, immune checkpoint blockade; IPA, ingenuity pathway analysis.
According to the workflow in Figure S1A, we identified the differentially expressed (DE) genes of R and NR in the 13 cell types using the scCODE R package (22) and obtained a total of 26 DE gene lists (Extended data 1, see Study Design for details). We sought to identify gene sets that were effective predictors of ICB outcomes. From the Molecular Signatures Database (MSigDB) (23), the 26 DE gene lists were found to be enriched in 1008 gene sets, including GOBP, Hallmark, KEGG, and Reactome (Extended data 2, see Study Design for details). The predictive capabilities of the 1008 gene sets were evaluated in terms of receiver operating characteristic (ROC) curves obtained with the Cancerclass R package (24) in the CheckMate cohort (20) (n.patients = 172, R = 39, NR = 133, see Methods for details). ROC p values for these 1008 gene sets were visualized and are listed in Figures S2 A–D and Extended data 2.
We found that the GOBP gene sets had overall better predictive performance than Hallmark, KEGG and Reactome. To identify the cell types with good predictive capability, we examined the enriched GOBP gene sets with ROC p values < 0.05 (without correction) within each DE gene list. A DE gene list was considered significant if no less than half of its enriched gene sets had ROC p values < 0.05. Next, based on the same criteria, significant DE gene lists were identified in Hallmark, KEGG, and Reactome. As a result, we identified 7 prediction-related DE gene lists (Figure S2A). The DE gene lists were further reduced to 3 with FDR corrected ROC p values (Figure S2E): the monocytes from Bi’s dataset responders (Bi’s_Mono.R), and CD4Ts and Tregs from Au’s dataset nonresponders (Au’s_CD4T.NR and Au’s_Treg.NR). We speculated that the CD4T, Treg and monocyte subtypes might be correlated with ICB responses.
MKI67+ CD4Ts and MKI67+ Tregs were enriched in nonresponders
The two significant nonresponder-associated DE gene lists: Au’s_CD4T.NR and Au’s_Treg.NR in Au’s dataset was further analyzed. The 1308 CD4Ts (NR.cells = 712, R.cells = 596) and 1501 Tregs (NR.cells = 126, R.cells = 375) were clustered into five CD4Ts subclusters and four Tregs subclusters by Seurat (25), respectively (Figures 2A, D). Subsequently, cluster-specific marker genes were identified by FindAllMarkers (Seurat), and the expression heatmap of the top 10 marker genes is shown in Figures 2B, E. The classical marker genes of CD4Ts and Tregs were highly expressed in all subclusters (Figures S3A, S3B).
Figure 2 Proliferative subtypes of CD4Ts and Tregs were enriched in nonresponders. Au’s CD4Ts and Tregs scRNA-seq datasets were analyzed. (A) UMAP plot of Au’s dataset’s CD4Ts that were classified into 5 subclusters from R and NR of distinct ICB outcomes. Bar plots show cell proportions grouped by clusters (left) and ICB outcomes (right). (B) Heatmap of scaled normalized expression for the top 10 specific marker genes of Au’s CD4T subclusters as identified by a two-sided Wilcoxon rank-sum test with FDR correction (q < 0.05). (C) ICB response prediction-related CD4T subcluster was identified by locating the effective (ROC p.adjust < 0.05) predictive gene set expression via gene set variation analysis (GSVA). (D) UMAP plot of Au’s Treg dataset that was classified into 4 subclusters of distinct immune checkpoint therapy outcomes from R and NR. Bar plots show cell proportions grouped by clusters (left) and ICB outcomes (right). (E) Heatmap of scaled normalized expression for the top 10 specific marker genes of Au’s Treg subclusters as identified by a two-sided Wilcoxon rank-sum test with FDR correction. (F) ICB response prediction-related Treg subtype was identified by locating the effective predictive gene set expression via GSVA. (G) Violin plot of the expression levels of proliferative marker genes in Au’s CD4Ts and Treg subclusters. A two-sided Wilcoxon test was used to determine significance between the subclusters of interest and others. ****P < 0.0001. (H) GOBP, Hallmark, KEGG and Reactome analysis results of Au’s CD4Ts subcluster 4 (CD4T_C4) and Au’s Treg subcluster 3 (Treg_C3) compared with other subclusters.
There were 7 gene sets of Au’s_CD4T.NR and 5 gene sets of Au’s_Treg.NR with ROC p.adjust < 0.05 (Figure S2E). We analyzed the expression of these gene sets via gene set variation analysis (GSVA) (26) and found that they were highly expressed in CD4T subcluster 4 (CD4T_C4, NR.cells = 24, R.cells = 8) (Figure 2C, Figure S3C) and Treg subcluster 3 (Treg_C3, NR.cells = 96, R.cells = 5) (Figure 2F, Figure S3D). Therefore, both CD4T_C4 and Treg_C3 cells were likely associated with ICB resistance.
Further analysis of these two subclusters revealed that although they belong to different cell types, they are both proliferative cell populations. Specifically, they both highly expressed MKI67, TOP2A, TUBB, TUBA1B, STMN1, TYMS and other proliferation markers (27, 28) (Figure 2G, Figures S3E, S3F). The GOBP, Hallmark, KEGG and Reactome analyses of their marker genes showed that cell cycle-related pathways and metabolic activities (e.g., DNA repair and oxidative phosphorylation [OXPHOS], etc.), oncogenic pathways such as MYC (29) and E2F (30) targets were dramatically enriched (Figure 2H, Extended data 3). In contrast, pathways associated with tumor suppression, such as immune activation-related pathways (31), antigen presentation and processing (16), IFNγ responses (15, 32), TNFα signaling (33), and P53 pathways (34), were remarkably suppressed in both subclusters (Figure 2H, Extended data 3). In addition, CD4T_C4 also downregulated anticancer immune responses and cytokine-related pathways (35), while Treg_C3 downregulated anticancer IFNα and IFNβ signaling (15, 32) and activated oncogenic mTORC1 signaling (36) (Figure 2H, Extended data 3). Their top 20 pathways sorted by the absolute value of the normalized enrichment score (NES) also confirmed the above observation (Figures S3G, S3H). Next, we used ingenuity pathway analysis (IPA) for further verification. IPA showed that CD4T_C4 and Treg_C3 cells have significantly activated cell cycle-related kinetochore, chromosomal replication, cyclins and mitosis, as well as oxidative phosphorylation, while cell cycle-related checkpoints were inhibited in CD4T_C4 cells (Figures S3I, S3J). Such observations were consistent with previous studies, i.e., dysfunction of cell cycle-related checkpoints often results in genomic instability and oncogenesis (37). Therefore, CD4T_C4 and Treg_C3 cells are proliferative subtypes and may impair antitumor immunotherapy. In view of the high expression of proliferation markers and enrichment of cell cycle-related pathways in both CD4T_C4 and Treg_C3, we annotated them as MKI67+ CD4Ts and MKI67+ Tregs, respectively.
Validation of MKI67+ CD4Ts and MKI67+ Treg signatures in independent datasets
As mentioned above, we preliminarily confirmed that MKI67+ CD4Ts and MKI67+ Tregs were related to ICB resistance. To further validate this finding, we investigated whether tumors enriched with these two cell subtypes were more susceptible to ICB resistance.
The 70 marker genes of MKI67+ CD4Ts (annotated as MKI67+CD4T.Sig) were characterized by cell cycle, MYC, E2F targets and OXPHOS (Figure S4A, Extended data 4). The GSVA (26) score of MKI67+CD4T.Sig was significantly higher than the other CD4Ts subclusters (Figure 3A, Figure S4B). These genes were enriched in nonresponders (Figure 3B), and the GSVA score was notably higher in nonresponders as well (Figure S4C). Subsequently, we validated MKI67+CD4T.Sig in the pretreatment CheckMate cohort (n.patients = 172, R = 39, NR = 133) (20) and found that nonresponders had higher GSVA scores (p = 0.028) (Figure 3C). Consistently, gene set enrichment analysis (GSEA) showed significant enrichment of MKI67+ CD4T.Sig in nonresponders compared to responders (Figure S4D).
Figure 3 Validation of the MKI67+ CD4Ts signature using independent bulk RNA-seq datasets. Au’s CD4Ts scRNA-seq dataset and CheckMate cohort were analyzed. (A) Violin plot of a 70-gene signature (MKI67+CD4T. Sig) GSVA scores showed that it can specifically characterize proliferative subcluster 4 of CD4Ts (MKI67+ CD4Ts). Centerline, median. Box limits, upper and lower quartiles. Whiskers, 1.5 interquartile range. Points beyond whiskers, outliers. A two-sided Wilcoxon test was used to determine significance. ****P < 0.0001. (B) GSEA revealed that MKI67+CD4T.Sig was enriched in NR of MKI67+ CD4Ts. The p value was FDR-adjusted by the FDR method. (C) Boxplot validated that MKI67+CD4T. The Sig GSVA scores of NR were significantly higher than those of R in the CheckMate cohort (n.patients = 172, R = 39, NR = 133) via GSVA analysis. Centerline, median. Box limits, upper and lower quartiles. Whiskers, 1.5 interquartile range. Points beyond whiskers, outliers. A two-sided Wilcoxon test was used to determine significance. (D), (E) Survival analysis was performed on GSVA scores for MKI67+CD4T cells. Sig in the pretreatment CheckMate cohort (n.patients = 181, R = 39, NR = 133, NE = 9, NEs were not removed for survival analysis) based on the Kaplan–Meier method. The groups were dichotomized at the median GSVA score, and the log-rank test was used to determine significance. Dashed line: median survival time. Color range: 95% confidence interval (CI).
Next, we analyzed the impact of MKI67+ CD4Ts on the prognosis of ICB treatment patients. It was found that in the pretreatment CheckMate cohort (n.patients = 181, R = 39, NR = 133, NE = 9, NEs were not removed for survival analysis, the same below), the group with low MKI67+CD4T.Sig GSVA scores were associated with significantly higher overall survival (OS, p = 0.022) and progression-free survival (PFS, p = 0.022) (Figures 3D, E). Subsequently, we corroborated this finding in the Javelin101 cohort (4) with 726 pretreatment bulk RNA-seq data (see Methods for details). Consistent results were obtained, i.e., the group with low GSVA scores had markedly longer PFS (p < 0.0001) (Figure S4E). To eliminate the effects of sunitinib treatment on the above results, we extracted and analyzed the avelumab plus axitinib group (4) bulk RNA-seq data (n.patients = 354), and the results were the same (p = 0.022, Figure S4F). In addition, we validated the results in the posttreatment bulk RNA-seq dataset of the CheckMate 009 (CM009) cohort (38) (CM009_POST, n.patients = 55, R = 5, NR = 47, NE = 3, see Methods for details). Lower MKI67+CD4T.Sig GSVA scores after ICB treatment were also significantly correlated with prolonged OS (p = 0.018) and PFS (p = 0.0093) (Figures S4G, S4H). Finally, univariate logistic regression with MKI67+CD4T.Sig in predicting ICB outcomes obtained an AUC of 0.62 (p = 0.008, Figure S4I).
Similar procedures were used to verify MKI67+ Tregs. The MKI67+Treg.Sig contained 80 genes involved in the cell cycle, MYC and E2F targets, mTORC1 signaling, OXPHOS and DNA repair (Figure S5A, Extended data 4). Violin and feature plots of GSVA scores showed that this gene set had distinguished specificity for characterizing MKI67+ Tregs (Figures S5B, S5C). In both single-cell and bulk RNA-seq datasets, nonresponders had higher MKI67+Treg.Sig GSVA scores than responders (p values 0.046 and 0.0047, respectively) (Figures S5D, S5E) and had significant MKI67+Treg.Sig enrichment (p < 0.001 and 0.0013, respectively) (Figures S5F, S5G). Survival analysis of the CheckMate, Javelin 101 and CM009_POST cohorts consistently demonstrated that high levels of MKI67+ Tregs, either pre- or post-ICB treatment, were all associated with worse OS and PFS (Figures S5H–S5L). The above results confirmed that MKI67+ CD4Ts and MKI67+ Tregs could shorten OS and PFS and promote ICB resistance. Similarly, univariate logistic regression showed that MKI67+Treg.Sig performed well in predicting ICB outcomes (AUC = 0.64, p = 0.011, Figure S5M).
Finally, we used an independent single-cell dataset, Bi’s dataset, to validate the universality of MKI67+ CD4Ts and MKI67+ Tregs characteristics. The 2245 CD4Ts (NR.cells = 818, R.cells = 641, NoICB.cells = 786) and 740 Tregs (NR.cells = 89, R.cells = 307, NoICB.cells = 344) from Bi’s dataset were clustered into five CD4Ts subclusters (Figure S6A) and three Tregs subclusters (Figure S6D), respectively. We matched subclusters with high levels of MKI67+ CD4Ts and MKI67+ Treg signatures through GSVA with CD4T_C0 (NR.cells = 56, R.cells = 347, NoICB.cells = 277) and Treg_C0 (NR.cells = 20, R.cells = 156, NoICB.cells = 106) (Figures S6B, S6E), respectively. As previously described, there were 7 gene sets of Au’s_CD4T.NR and 5 gene sets of Au’s_Treg.NR with ROC p.adjust < 0.05 (Figure S2E). In parallel, we found that these gene sets were also specifically enriched in CD4T_C0 and Treg_C0 of Bi’s dataset (Figures S6C, S6F), respectively. Such consistency was also confirmed by overall expression (OE) (39). The results showed that nonresponders in CD4T_C0 had higher MKI67+ CD4Ts-OE signals than responders (p = 0.018, Figure S6G), and a similar result was observed in Treg_C0 cells (p = 0.0004, Figure S6H).
Antigen-presenting monocytes are associated with the ICB response
As mentioned before, some gene sets enriched in responder monocytes (Bi’s_Mono.R) had good predictive capability (Figure S2E, Extended data 2), which suggested that certain subtypes of monocytes might be associated with the ICB response. The UMAP of monocytes (NR.cells = 122, R.cells = 364, NoICB.cells = 669) in Bi’s dataset showed four subclusters (Figures 4A, B). The canonical marker genes of monocytes were highly expressed in all subclusters (Figure S7A).
Figure 4 Antigen-presenting monocyte subtype promotes ICB response and is validated in independent bulk RNA-seq datasets. Bi’s monocyte scRNA-seq dataset and the CheckMate cohort were analyzed. (A) UMAP plot of monocytes from Bi’s dataset that were classified into 4 subclusters from ICB-R, ICB-NR and NoICB samples. Bar plots show cell proportions grouped by clusters (left) and ICB outcomes (right). (B) Heatmap of scaled normalized expression for the top 10 specific marker genes of Bi monocyte subclusters as identified by a two-sided Wilcoxon rank-sum test with FDR correction (q < 0.05). (C) The ICB response prediction-related monocyte subcluster was identified by locating the effective predictive gene set expression via GSVA. (D) GOBP, Hallmark, KEGG and Reactome analysis results of Bi’s monocyte subcluster 0 compared with other subclusters. The top 20 pathways sorted by the absolute value of the normalized enrichment score (NES) are listed. (E) Violin plot of the expression levels of MHC class II genes in Bi’s monocyte subclusters. A two-sided Wilcoxon test was used to determine significance between the subclusters of interest and others. ****P < 0.0001. (F) Feature plot of a 45-gene signature (Mono_C0. Sig) GSVA scores showed that it can specifically characterize the antigen-presenting subcluster 0 of monocytes (Mono_C0). (G) GSEA showed that Mono_C0. Sig was enriched in R of Mono_C0. The p value was FDR-adjusted by the Benjamini–Hochberg method. (H) Boxplot validated that Mono_C0. The Sig GSVA scores of R were significantly higher than those of NR in the CheckMate cohort (n.patients = 172, R = 39, NR = 133) via GSVA analysis. Centerline, median. Box limits, upper and lower quartiles. Whiskers, 1.5 interquartile range. Points beyond whiskers, outliers. A two-sided Wilcoxon test was used to determine significance. (I) Kaplan–Meier plot of progression-free survival (PFS) data for the pretreatment CheckMate cohort (n.patients = 181, R = 39, NR = 133, NE = 9) on the basis of Mono_C0. Sig GSVA scores. The groups were dichotomized at the median GSVA score, and the log-rank test was used to determine significance. Dashed line: median survival time. Color range: 95% CI.
There were five gene sets with ROC p.adjust < 0.05 in Bi’s_Mono.R (Figure S2E), and they were uniquely highly expressed in subcluster 0 (Mono_C0, NR.cells = 56, R.cells = 262, NoICB.cells = 134) (Figure 4C, Figure S7B). Mono_C0 highly expressed MHC class II antigen and IFNγ response-related genes, including HLA-DRB1, HLA-DPA1, HLA-DPB1, CD74 and TUBA1B (Figures 4B, E). GSEA showed that Mono_C0 promoted the ICB response through multiple pathways, including antigen processing and presentation (16), adaptive immune response (31) and IFNγ signaling pathways (40) (Figure 4D, Extended data 3). In addition, IPA showed significant inhibition of the PD-1/PD-L1 pathway (Figure S7C), indicating that Mono_C0 may play an important role in immunotherapy.
To validate the role of antigen-presenting monocytes in immunotherapy, a Mono_C0.Sig gene panel was constructed based on its marker genes. This gene panel consisted of 45 genes, all of which represented characteristic pathways such as antigen processing and presentation, IFNγ signaling, and immune response (Figure S7D, Extended data 4), with good specificity for characterizing Mono_C0 (Figure 4F, Figure S7E) and significant enrichment in responders (Figure 4G, Figure S7F). We then validated Mono_C0.Sig in the CheckMate cohort (n.patients = 172). GSVA and GSEA analyses revealed significant enrichment in responders (Figure 4H, Figure S7G). Survival analysis showed that higher GSVA scores indicated better PFS (p = 0.042) (Figure 4I), although the GSVA scores were not correlated with OS (p = 0.95) (Figure S7H). Consistently, the same results were observed in the CM009_POST cohort (Figures S7I, S7J). The results suggested that ccRCC patients with abundant Mono_C0 subtypes had prolonged PFS after ICB treatment. Such observations were consistent with previous reports, e.g., MHC class II molecule expression was related to the ICB response and prognosis improvement (41, 42). Furthermore, univariate logistic regression analysis showed that Mono_C0.Sig could predict ICB outcomes with an AUC of 0.63 (p = 0.019, Figure S7K).
Comparing the NoICB and ICB treatment samples from Bi’s dataset, the expression of Mono_C0. Sig was elevated after ICB treatment (p < 0.001) (Figure S7L). We further analyzed a subdataset from the CM009 cohort (CM009_Paired), which contained paired kidney biopsy samples of ICB pretreatment and week 4 treatment from 42 patients (38). A similar dynamic pattern was observed (Wilcoxon’s paired test p = 0.0068) (Figure S7M). Evaluation in the CM009_POST cohort indicated the upregulation of Mono_C0.Sig was more obvious in ICB responders than nonresponders (p = 0.038) (Figure S7N).
ICB response prediction model construction and validation
As MKI67+ CD4Ts and MKI67+ Tregs tended to impair immunotherapy efficacy and antigen-presenting monocytes promoted the ICB response, we compared the predictive capability of these three cell subtypes with that of other CD4Ts, Tregs and monocytes subclusters. Specifically, the top 500 genes sorted in ascending order of adjusted p value (obtained by FindAllMarkers of Seurat) were selected to identify the top 10 enriched GOBP gene sets for each cell subcluster of CD4Ts, Tregs and monocytes (Extended data 5, Figure S8). We subsequently tested the predictive capability of each individual gene set in the CheckMate cohort (n.patients = 172) using Cancerclass. Most gene sets associated with MKI67+ CD4Ts, MKI67+ Tregs and antigen-presenting monocytes had high predictive capability compared with the rest (Figure S9A-S9C, Extended data 5). This provided a foundation to develop a signature panel for ICB outcome prediction based on the three cell subtypes.
Based on the selected gene sets with significant predictive capability (p.adjust < 0.05) (Figures S9A-S9C), we developed the overall gene signatures using Cancerclass according to the workflow shown in Figure S8 (see Methods for details). Initially, 209 genes were obtained (Figure S8, orange box), and GSEA showed that the top 20 pathways were mainly the cell cycle, DNA metabolism and repair, E2F and MYC targets (Figure S10A). Most of these pathways contributed to the poor immune response according to previous studies (29, 30, 43). The top 10 pathways were used to predict the outcomes of the CheckMate cohort, and these pathways had good prediction performance, with AUCs between 0.68 and 0.83 (Figure S10B). To construct a more effective prediction model, we executed the cycle algorithm shown in Figure S8 (purple box, see Methods for details) and examined the AUCs for combinations with different numbers of genes (Figure 5A, Extended data 6). A good prediction model should balance the highest reliability and fewest genes. Therefore, to balance these two factors, we selected the 47-gene combination as our ccRCC ICB response prediction signature - ccRCC.Sig in downstream analysis (Figure 5A, dotted line). Details of these 47 genes can be found in Supplementary Table 2. The ccRCC. Sig could accurately discriminate between responders and nonresponders in the CheckMate cohort (n.patients = 172, R = 39, NR = 133). Specifically, the AUC of this signature was 0.93 (95% confidence interval [CI]: 0.91-0.95), the sensitivity was 80% (74-86%), the specificity was 92% (81-98%), and p = 7.2e-07 (Figure 5B). This response prediction signature gathers characteristic genes from MKI67+ CD4Ts, MKI67+ Treg cells and antigen-presenting monocyte subtypes (Figures S11A–C). In addition, we tested the performance of 30-gene and 20-gene combinations in the CheckMate cohort. However, their predictive performances were less than that of ccRCC.Sig, they were also effective, with AUCs of 0.89 (0.87-0.91) and 0.86 (0.84-0.89), respectively (Figures S12A, B).
Figure 5 The ccRCC.Sig signature can effectively predict ICB outcomes in ccRCC patients. The bulk RNA-seq datasets CheckMate (n = 172, R = 39, NR = 133), CM025 (n = 111, R = 25, NR = 86), CM010 (n = 45, R = 11, NR = 34), CM009_PRE (n = 56, R = 9, NR = 47), PMC6035749 (n = 17, R = 5, NR = 12) and GSE67501 (n = 11, R = 4, NR = 7) were analyzed. (A) Bar graph showing the AUC of gene combinations with maximal AUC for each cycle (different gene-number combinations). Dotted line: 47-gene combination - ccRCC.Sig. (B) ccRCC.Sig had significantly high predictive value for ICB outcomes in the CheckMate cohort. (C), (D) ccRCC.Sig accurately predicted the ICB outcomes in the validation (C) PMC6035749 and (D) GSE67501 cohorts. (E) Comparison of the performance (AUC and p value) of ccRCC.Sig with 13 other ICB response signatures. (F) Comparison of the performance (AUC) of ccRCC.Sig with other clinical traits and molecular features in the CheckMate cohort (115 samples remained after deleting missing data samples). (G) Validation of ccRCC.Sig using three other machine learning algorithms. SVM, support vector machine; NB, naïve Bayes; KNN, k-nearest neighbors. (H) Univariate Cox regression analyses of ccRCC.Sig, clinical traits and molecular features. TMB, tumor mutational burden; ITH, intratumor heterogeneity.
Subsequently, we validated the performance of ccRCC.Sig in independent pretreatment bulk RNA-seq datasets. Because the CheckMate cohort (20) is integrated data of the CheckMate 025 (CM025) (5), CheckMate 010 (CM010) (44), and CM009 (45) cohorts (see Methods for details), we first validated separately in each of these three cohorts. For the CM025 (n.patients = 111, R = 25, NR = 86) and CM010 cohorts (n.patients = 45, R = 11, NR = 34), the AUCs of this signature were 0.88 (0.85-0.91) and 0.99 (0.98-1), respectively (Figures S12C, S12D). For the CM009 cohort, we used its larger pretreatment dataset (38) (CM009_PRE, n.patients = 56, R = 9, NR = 47), and this signature possessed an AUC of 0.89 (0.85-0.93) (Figure S12E). In addition, two more independent ccRCC bulk RNA-seq datasets were utilized: PMC6035749 (46) and GSE67501 (47). For PMC6035749 (n.patients = 17, R = 5, NR = 12), the performance of ccRCC.Sig was AUC = 0.87 (0.79-0.96), sensitivity = 100% (78-100%), and specificity = 60% (19-92%) (Figure 5C). For GSE67501 (n.patients = 11, R = 4, NR = 7), ccRCC.Sig accurately distinguished all responders from nonresponders, with AUC = 1 (1–1), sensitivity = 100% (65-100%), and specificity = 100% (47-100%) (Figure 5D).
To further justify the predictive capability of ccRCC.Sig, we compared ccRCC.Sig with 13 other ICB response signatures previously reported (8, 48–58), including the recognized IMPRES signature (51). The details of these signatures are listed in Supplementary Table 3. They were used as independent classifiers to predict the ICB outcomes of six bulk RNA-seq cohorts using Cancerclass. The results show that ccRCC.Sig dramatically improved the performance of ICB response prediction in most cohorts (Figure 5E, Figures S12F-S12I, Supplementary Table 4). As a comparison, the performance of IMPRES.Sig ranked fourth in both the CheckMate and CM009_PRE cohorts (Figure 5E, Figure S12I) and fifth and third in the CM025 and CM010 cohorts (Figures S12G, S12H). We also compared ccRCC.Sig with clinical traits (e.g., age) and other molecular features (e.g., TMB, PD-1, tumor purity and angiogenesis, etc.), and the results showed the AUC of ccRCC.Sig was 0.91 (115 samples remained after deleting missing data samples) and that of other features was 0.48–0.62 (Figure 5F). In addition, we validated ccRCC.Sig using three more machine learning algorithms, i.e., support vector machine (SVM), naïve Bayes (NB) and k-nearest neighbors (KNN). As shown in Figure 5G, ccRCC.Sig still performed well in most cohorts. These results indicate that its predictive performance is stable.
In particular, prediction using Cancerclass is able to generate continuous prediction scores (z score) based on gene expression levels and convert them to probabilities (24). Therefore, to assess a patient’s resistance risk after ICB treatment, we converted the prediction scores to nonresponse probabilities and estimated the risk by logistic regression (Figure S12J). Based on the pretreatment RNA-seq data of tumor patients, we can estimate the resistance probability and provide guidance and a reference for patients to decide whether to accept ICB treatment by using this model.
Moreover, we explored the relationship between ccRCC.Sig and the prognosis of patients. Univariate Cox regression demonstrated the pretreatment GSVA score of ccRCC.Sig, PBRM1 mutation and intratumor heterogeneity (ITH) were closely related to poorer OS (Figure 5H), and multivariate Cox regression found that ccRCC.Sig, wild-type PBRM1 (20, 46, 59) and ITH were independent risk factors (Figure S12K).
We also found a relationship between ccRCC.Sig and a few known ICB response factors. The results showed that ccRCC.Sig was positively correlated with immune checkpoint molecules (PD-1, CTLA-4 and IDO1, Figures S13A–S13C) but not with TMB or ITH (Figures S13D, S13E). This is consistent with previous studies. Although TMB and ITH are associated with the ICB response in various tumors (11, 60–62), they are not associated with the ICB response in ccRCC (4, 7, 12, 13, 20). Angiogenesis (63, 64) and myeloid infiltration (7, 63) have been reported to be related to the ICB response and resistance in ccRCC, respectively. Our results show that ccRCC.Sig was negatively correlated with pro-response angiogenesis signaling (Figure S13F) but positively correlated with pro-resistant myeloid infiltration (Figure S13G).
Discussion
Herein, we reanalyzed two publicly available ccRCC scRNA-seq datasets (12, 14) to explore effective predictive immunocyte subtypes and signatures in ccRCC. We identified three cell subtypes that were closely related to ICB outcomes: MKI67+ CD4Ts and MKI67+ Tregs contribute to ICB resistance, while antigen-presenting monocytes are correlated with the ICB response. We harnessed the gene markers associated with the above three subtypes and developed a signature, ccRCC.Sig. The prediction capability of ccRCC.Sig was systematically evaluated using various datasets, modeling methods and known risk factors. Our analyses demonstrated that compared to conventional prediction factors, ccRCC.Sig dramatically improved the reliability of ICB outcome prediction for ccRCC patients.
There are few ICB response prediction biomarkers currently available for ccRCC, and the predictive performance of these biomarkers varies among different tumor types (65). Immunohistochemistry (IHC) to detect PD-L1 expression on tumor cells or tumor-infiltrating immunocytes is the first clinically validated and the most widely used biomarker currently in ICB therapy (66). However, a lack of PD-L1 expression is considered an insufficient negative predictor of immunotherapy response, and patients with IHC PD-L1 (–) may also benefit from ICB in some clinical trials (67). Another biomarker is TMB, which is a measure of total tumor mutation. However, the clinical application of TMB detection is limited by high technical requirements, complex data and the need for bioinformatics experts. Next-generation sequencing (NGS) is immature and expensive. In addition, flow cytometry can be used for the development and detection of biomarkers. However, this method also has limitations of high cost, complex operations, high technical requirements and limited antibody availability in clinical application. The predictive biomarker we developed, ccRCC.Sig has high sensitivity and specificity and can be detected by real-time quantitative PCR, which enables it to be a highly reliable and clinically practical prediction tool.
Studies have pointed out that CD4+ T cells are necessary for effective antitumor immunity (68). CD4Ts exert antitumor activity through various direct and indirect mechanisms, including cytolysis of tumors and modulation of the tumor microenvironment (TME) (69). It can increase the activity and quality of B cells and cytotoxic T lymphocyte (CTL) responses via cell–cell and cell surface receptor interactions (70). Tregs, although a major subset of CD4Ts, differ from traditional CD4Ts and normally mediate immunosuppression and tolerance in homeostasis and inflammation (71). The consensus holds that Tregs in the TME exert an inhibitory effect against tumor immunity (72). In this study, MKI67+ CD4Ts and MKI67+ Tregs were identified as belonging to proliferative subgroups. They are two cell subtypes with similar specific functional characteristics, which are also similar to the proliferative CD8Ts described by Borcherding et al. (27). Although the effects of proliferative MKI67+ CD4Ts/Tregs on ICB outcomes have been reported in other tumor types (73, 74), this correlation has not been documented in ccRCC before.
In recent years, researchers have identified similar proliferative populations in many tumor types, such as melanoma (28, 48), lung cancer (75, 76), gastric cancer (73), head and neck cancer (74), and colon cancer (77), by scRNA-seq analysis or other experimental methods. The characteristics of these populations appear to be shared across different tumors (28), including ccRCC (27, 78). However, studies targeting these proliferative populations have mostly focused on CD8Ts. Interestingly, a proportion of hyperproliferative cells were also present in Tregs and CD4Ts, similar to exhausted/dysfunctional CD8Ts (28, 75). These CD8Ts have several common features, including active proliferation (MKI67+), high immune checkpoint or coinhibitory receptor expression (PD-1+, LAG3+, CTLA4+), terminal differentiation (EOMES+) and defective IFNγ production (28, 75, 76, 78). In esophageal squamous cell carcinoma, the proliferative CD4-C5-STMN1 cell subtype has the most exhausted features (79). Studies have shown that PD-1+ CD8Ts proliferate explosively after PD-1 blockade and impair the ICB response in chronic viral infections and cancers (76, 80). Similarly, a study of sarcoidosis showed that PD-1 blockade rescued the proliferation of CD4Ts (81). We also examined the expression of the above proliferative CD8Ts-associated characteristic genes in our MKI67+ CD4Ts and MKI67+ Tregs. We found that in these two cell subtypes, some immune checkpoint molecules, effector molecules and exhausted marker genes were preferentially enriched (Figures S14A, S14C). This finding indicates that they also appear to share these dysfunctional features with proliferative PD-1+ CD8Ts.
We also found that MKI67+ CD4Ts and MKI67+ Tregs belong to the PD-1+ cell population (Figures S14B, S14D). Consistently, a portion of Tregs have recently been found to express moderate levels of PD-1 (73, 74), in addition to activated and exhausted CD8Ts. PD-1+ Tregs show high proliferation and DNA metabolism features in head and neck cancer (74), and the proliferative activity of these cells with hyperprogressive disease (HPD) is significantly increased after ICB treatment in gastric cancer and compromises ICB treatment (73). PD-1 expressed in these cells might be a negative regulator of proliferation and immunosuppressive activity, while PD-1 blockade enhanced those of these cells.
Monocytes bridge the innate and adaptive immune response, with a paradoxical effect on tumors, involving multiple mechanisms, such as immune tolerance, angiogenesis, tumor cell metastasis, antitumor effector production, and antigen-presenting cell activation (82). These effects depend on monocyte plasticity under TME stimulation. It has been reported that circulating monocytes can infiltrate mucosa or inflammatory tissues and differentiate into monocyte-derived macrophages (mo-Macs) or monocyte-derived dendritic cells (mo-DCs) (83). Studies have shown that mo-Macs are powerful producers of proinflammatory cytokines (e.g., TNFα, interleukin, and IFNγ) (84), which are essential for immunocyte recruitment and immune response triggering. In addition, mo-Macs isolated from ascites are capable of cross-presenting antigens via specific pathways (85). In addition to MHC class II molecular genes, the expression levels of proinflammatory cytokines such as TNF, IFNγ, and IL18, as well as complement system genes such as C1QA, C1QB, and C3, were also significantly higher in our antigen-presenting monocytes (Mono_C0) than in the other monocyte subclusters (Figure S14E). Antigen processing and presentation, IFNγ signaling, and inflammatory response/regulation-related pathways were also dramatically enriched in this cell subtype (Figure 4D, Extended data 3). Based on the above characteristics, we speculated that Mono_C0 probably belongs to a monocyte subgroup similar to functional mo-Macs. This hypothesis is also supported by Lam et al. (86), who found that microbial signals program monocytes to inflammatory macrophages, characterized by high expression of MHC class II molecules and intermediate levels of F4/80 and CD68. This microbial-induced transformation can effectively improve the ICB response rate of melanoma.
In summary, our study not only provides an effective prediction model for ccRCC immunotherapy but also provides a pipeline for the development of therapy prediction models. Our method is applicable for prediction model development for various tumor therapies. There are still limitations in this study. We only described the relationship between the above three cell types and ICB outcomes but did not clarify their mechanism. In the future, large sample experiments with more rigorous designs are needed to explore the mechanisms of these cells to consolidate the biological findings in this study.
Methods
scRNA-seq and bulk RNA-seq datasets collection
To explore immunocyte subtypes and signatures with good predictive capability for ICB outcomes, we downloaded and reanalyzed two available ccRCC scRNA-seq datasets published in the initial article (12, 14) and validated the results in multiple ccRCC pretreatment bulk RNA-seq datasets. The first scRNA-seq dataset (Bi’s dataset) (14) was obtained from fresh biopsy or surgical resection samples of 8 patients, of which 5 patients were treated with ICB (anti-PD-1 combined with or without TKIs/anti-CTLA-4), and the other 3 patients did not receive systemic treatment. According to RECIST v1.1 (21), two of the five ICB-treated patients had PR efficacy, and the other three had SD, PD, and NE efficacy. In this study, the NE sample (P912) was omitted throughout. Bi’s dataset is available from the Single Cell Portal (https://singlecell.broadinstitute.org/single_cell/study/SCP1288/tumor-and-immune-reprogramming-during-immunotherapy-in-advanced-renal-cell-carcinoma#study-summary). The second scRNA-seq dataset (Au’s dataset) (12) was obtained from nephrectomy samples from 2 patients treated with nivolumab with RECIST efficacy of PD and PR. The single-cell count matrices and metadata can be downloaded publicly from the University College London (UCL) website (https://doi.org/10.5522/04/16573640.v1). Detailed clinicopathological information for all patients in Bi’s dataset and Au’s dataset, as well as the processing, clustering and cell type definition methods of the scRNA-seq datasets, were described in detail in their original articles (12, 14).
The ICB pretreatment bulk RNA-seq data used for validation consisted primarily of CheckMate 009 (CM009, NCT01358721) (45), CheckMate 010 (CM010, NCT01354431) (44), CheckMate 025 (CM025, NCT01668784) (5) and Javelin101 (NCT02684006) (4) cohorts. The first three are prospective phase I, phase II, and phase III clinical trials of nivolumab in the treatment of advanced ccRCC. The CheckMate cohort was an integrated dataset of the above three bulk RNA-seq datasets selected against certain criteria, leaving 181 nivolumab-treated samples after the removal of the everolimus-treated samples. This cohort included 16 CM009 samples (R = 3, NR = 13), 45 CM010 samples (R = 11, NR = 34), and 120 CM025 samples (R = 25, NR = 86, NE = 9). The selection strategy and the integration, alignment, quantification, and batch effect correction methods for this dataset have been described in detail by Braun et al. (20), and the normalized RNA-Seq expression data are provided in their supplementary information. The Javelin101 cohort is a randomized phase III clinical trial (NCT02684006) of avelumab (anti-PD-L1) plus axitinib (TKI) versus sunitinib (multitarget TKI) for the treatment of ccRCC (4). This cohort includes pretreatment bulk RNA-seq data from a total of 726 patients, 354 of whom received avelumab plus axitinib and 372 of whom received sunitinib, and its normalized RNA-seq data can be found in the supplementary original article. In addition, we downloaded the CM009 cohort dataset, which has a larger number of patients (38) and can be subdivided into three datasets: a pretreatment dataset (CM009_PRE, n.patients = 59, R = 9, NR = 47, NE = 3), a week 4 treatment dataset (CM009_POST, n.patients = 55, R = 5, NR = 47, NE = 3), and a pretreatment and week 4 treatment paired samples dataset (CM009_Paired, n.patients = 42, R = 5, NR = 37). They are publicly available from ArrayExpress (accession number: E-MTAB-3218). Two other independent datasets of anti-PD-1 treatment were also downloaded to test the predictive performance of ccRCC.Sig, GSE67501 cohort (47), available on Gene Expression Omnibus (GEO), and PMC6035749 cohort (46), available from supplementary of original article. The database used in our study included MSigDB v7.5.1 (23) (http://www.gsea-msigdb.org/gsea/index.jsp), from which the Hallmark, KEGG, GOBP, and Reactome gene sets were downloaded for GSEA.
Detailed metadata for all scRNA-seq and bulk RNA-seq dataset samples used in this study can be found in Supplementary Table 1.
Study design
According to the flow shown in Figure S1A, we extracted a total of 13 cell types from two scRNA-seq datasets. We then used scCODE v1.0.0.0 R package provided by Zou et al. (22) to identify the DE genes of responders and nonresponders in the above 13 cell types, respectively. scCODE can check the selected DE genes through a variety of testing methods, improving the accuracy of single-cell DE analysis. In this way, two DE gene lists, which were remarkably highly expressed in R and NR, were obtained for each cell type, with a total of 26 DE gene lists (Extended data 1. We then used the Investigate Gene Sets tool (http://www.gsea-msigdb.org/gsea/msigdb/annotate.jsp) to compute enriched gene sets between our gene lists and gene sets in MSigDB (23). Due to this tool’s limitation on the submitted gene numbers (≤ 500), for the 26 DE gene lists identified by scCODE, we sorted them according to the absolute value of logFC of each DE gene list from large to small. The gene sets enriched in GOBP, Hallmark, KEGG and Reactome were identified by submitting the top 500 genes (all genes were submitted if less than 500). With the default settings of this tool (show top 10 genes & FDR q-value less than 0.05), each DE gene list can obtain several gene sets enriched in GOBP, Hallmark, KEGG and Reactome. The predictive capability of these gene sets was tested using Cancerclass v1.34.0 R package (24). Due to the limitation of this package on the input gene numbers (minimum 3 genes), we excluded gene sets with fewer than 3 genes. Eventually, 1008 gene sets were obtained for subsequent analysis (Extended data 2).
We tested the predictive capability of the above 1008 gene sets on the ICB outcomes of the CheckMate cohort (n.patients = 72, R = 39, NR = 133) using the Cancerclass R package. The prediction sensitivity and specificity were assessed by ROC curve and corresponding AUC. This R package is dedicated to developing and validating classification tests for high-dimensional molecular data. Feature selection and nearest centroid classification were performed in sequence. Finally, the classification results were carefully verified by using the multiple random validation protocol to generate continuous prediction scores (24). Each gene set was tested as an independent classifier, and the p value of the ROC curve was calculated by Welch’s t test built into the Cancerclass R package, which reflects the effectiveness of the classifier’s classification test results. These 1008 ROC p values (Extended data 2) were used for subsequent sorting of cell subtypes.
Single-cell RNA sequencing data processing
We directly extracted the cell types defined by the authors in the original articles of Bi’s dataset (14) and Au’s dataset (12). For subtype analysis of each cell type, we performed fine clustering utilizing Seurat v4.0.4 R package (25). Specifically, we first performed across-sample integration on the extracted datasets. Prior to this, the expression of each gene was normalized by the total expression in the corresponding cells, multiplied by a scaling factor of 10,000, and then log2 transformed (Seurat default method). Two thousand variable features were identified by FindVariableFeatures and used for subsequent analysis. Then, the mutual nearest neighbors (MNN) method of the batchelor v1.0.1R package (87) was used to correct batch effects, and ScaleData were used to regress the percentage of mitochondrial transcripts. The integrated assay was used only for dimension reduction and clustering, and the raw log-normalized expression data were used for all DE and gene level analyses. Then, principal component analysis (PCA) was performed on the integrated assay by using RunPCA. The first 20 principal components were taken for Louvain clustering of cells with a resolution parameter of 0.5. Finally, visualization was performed in two-dimensional space by uniform manifold approximation and projection (UMAP) (Dims = 1:20). This procedure was used for scRNA-seq datasets from all cell type analyses.
Gene differential expression analysis
Using the FindAllMarkers built in the Seurat package, DE analysis was performed on each Louvain cluster and all other clusters by setting the parameters min.pct = 0.1 and logfc.threshold = 0.25. Genes with p.adjust < 0.05 were selected as cluster-specific marker genes. Using the FindMarkers built into the Seurat package, we compared all gene expression fold changes in nonresponders and responders in corresponding cells by setting the parameter logfc.threshold = -Inf, min.pct = -Inf, min.diff.pct = -Inf. Utilize the default settings of scCODE v1.0.0.0 The R package (22) (https://github.com/XZouProjects/scCODE) was used to analyze the DE genes of responders and nonresponders in each cell type. The generalized linear model method of Limma v3.46.0 The R package (88) (https://doi.org/doi:10.18129/B9.bioc.limma) was used to compute the gene expression fold changes of nonresponders and responders in the CheckMate cohort. We fit the linear model of each gene using the weighted least square method through the lmFit function and compared each gene through the contrasts.fit function. Finally, we used the empirical Bayes smoothing of standard errors. In all DE analyses, genes with Bonferroni FDR-corrected p values < 0.05 were considered DE genes by the bilateral Wilcoxon rank sum test.
GSVA, GSEA and IPA
We applied the GSVA method with default settings to assign a specific gene signature activity score for individual cells or samples, as implemented in GSVA v1.38.2 R Package (26) (https://doi.org/doi:10.18129/B9.bioc.GSVA). GSEA was performed on preordered DE gene lists based on predownloaded Hallmark, KEGG, GOBP, and Reactome gene sets using the default parameters of clusterProfiler v3.18.1 R package (89) (https://doi.org/doi:10.18129/B9.bioc.clusterProfiler). This R package can also be used to examine whether a particular gene set is enriched at the top or bottom of a preordered gene list. Gene sets with FDR-corrected p values < 0.05 by the Benjamini–Hochberg method were considered to be significantly enriched in one group when two groups were compared. Ingenuity pathway analysis was performed using QIAGEN IPA software (IPA Spring Release April 2022, [https://digitalinsights.qiagen.com/products-overview/discovery-insights-portfolio/analysis-and-visualization/qiagen-ipa/]) with lists of specific marker genes in the corresponding clusters as detected by Seurat.
Survival analysis, Cox and logistic regression analysis
Kaplan–Meier method-based survival analysis and univariate and multivariate Cox regression were performed on GSVA scores for specific gene signatures using Survival v3.2-13 R package and SurvivMiner v0.4.9 R package. To define “high-score” and “low-score” groups, in the CheckMate cohort, the groups were dichotomized at the median GSVA score. For other bulk RNA-seq datasets, we calculated the optimal cutpoint using the surv_cutpoint function built in the Survminer R package. The groups were compared using the log-rank test. Logistic regression analysis was performed using the R software embedded glm function for ICB outcomes and mean expression of specific gene signatures. The corresponding ROC curves were drawn by the pROC v1.18.0 R package.
Workflow of ICB response prediction signature development
Based on the cluster-specific genes of MKI67+ CD4T, MKI67+ Treg and Mono_C0 subtypes, we developed our gene signatures using Cancerclass v1.34.0 R package according to the flow shown in Figure S8. Specifically, the p values of ROC curves (30 in total) using the enriched GOBP gene sets of these three cells were obtained using Cancerclass v1.34.0 R package and then FDR-adjusted with the Benjamini–Hochberg method (Extended data 5). We united all the gene sets with p.adjust < 0.05 (26 gene sets in total) to obtain a total of 493 genes (union.genelist). Meanwhile, we compared all gene expression fold changes in NR and R through DE analysis in the corresponding cell subtypes. For MKI67+ CD4Ts, we selected the NR vs. R genes avg_log2FC > -0.1; for MKI67+ Tregs, we selected the NR vs. R genes avg_log2FC > -0.3; and for Mono_C0, we selected the NR vs. R genes avg_log2FC < 0.1. After intersecting these selected genes, we obtained a gene list of 11,884 genes (intersect.genelist). Subsequently, we took the intersection of union.genelist and intersect.genelist and excluded absent genes in the CheckMate cohort. We finally obtained the “genelist”, which contained 209 genes. Next, the genelist was used to execute the cycle algorithm shown in Figure S8 (purple box). First, we randomly selected combinations containing 208 genes from the genelist, resulting in 209 different combinations. The predictive capability of all these combinations was then tested, and the AUC was estimated in the CheckMate cohort using the Cancerclass R package. The gene combination with the largest AUC (genelistmaxAUC) of these 209 combinations was preserved and used for the next cycle. The cycle was repeated until the total gene number of the gene combination was < 3, and finally, all the circulating genelistmaxAUCs and their corresponding AUCs were output (Extended data 6). A combination with appropriate gene numbers was selected as ccRCC.Sig for subsequent analysis.
Statistical analysis
The predictive capability of each classifier in this study for ICB outcomes was assessed by plotting the ROC, calculating the AUC, and estimating the sensitivity and specificity implemented in Cancerclass v1.34.0 R package24, and its effectiveness was assessed by the ROC curve’s p value. The 95% confidence intervals for sensitivity and specificity were calculated by Wilson’s method built into Cancerclass v1.34.0 R package, and p-values were calculated by Welch’s t-test. Unless otherwise noted, all p-values in this manuscript were adjusted by the Benjamini–Hochberg method, and adjusted p-values < 0.05 were considered statistically significant. Wilcoxon’s test was used to classify variables. All confidence intervals were reported as binomial 95% confidence intervals. All statistical analyses in this study were performed with R v3.5.3 software.
Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors. The computer R codes for processing and analysis of this study are available at https://github.com/XZouProjects/ccRCC.Sig.
Author contributions
Conceptualization: XH and XZ. Methodology: XZ and JH. Investigation: XZ and KZ. Project administration: XH and XZ. Writing – original draft: KZ, LG, JH, XZ and XH. All authors contributed to the article and approved the submitted version.
Funding
This work was supported in part by the National Natural Science Foundation of China [82170045 to JH]; the Innovative Research Team of High-level Local Universities in Shanghai [SHSMU-ZLCX20212301 to JH] and the Shanghai Municipal Health Commission Project [202040375 to XH].
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2022.972227/full#supplementary-material
References
1. Saad A, Gad M, Al-Husseini M, Ruhban I, Sonbol M, Ho T. Trends in renal-cell carcinoma incidence and mortality in the united states in the last 2 decades: A seer-based study. Clin Genitourin Cancer (2019) 17(1):46–57.e5. doi: 10.1016/j.clgc.2018.10.002
2. Smittenaar C, Petersen K, Stewart K, Moitt N. Cancer incidence and mortality projections in the uk until 2035. Br J Cancer (2016) 115(9):1147–55. doi: 10.1038/bjc.2016.304
3. Şenbabaoğlu Y, Gejman R, Winer A, Liu M, Van Allen E, de Velasco G, et al. Tumor immune microenvironment characterization in clear cell renal cell carcinoma identifies prognostic and immunotherapeutically relevant messenger RNA signatures. Genome Biol (2016) 17(1):231. doi: 10.1186/s13059-016-1092-z
4. Motzer RJ, Robbins PB, Powles T, Albiges L, Haanen JB, Larkin J, et al. Avelumab plus axitinib versus sunitinib in advanced renal cell carcinoma: Biomarker analysis of the phase 3 javelin renal 101 trial. Nat Med (2020) 26(11):1733–41. doi: 10.1038/s41591-020-1044-8
5. Motzer R, Tannir N, McDermott D, Arén Frontera O, Melichar B, Choueiri T, et al. Nivolumab plus ipilimumab versus sunitinib in advanced renal-cell carcinoma. N Engl J Med (2018) 378(14):1277–90. doi: 10.1056/NEJMoa1712126
6. Giraldo N, Becht E, Pagès F, Skliris G, Verkarre V, Vano Y, et al. Orchestration and prognostic significance of immune checkpoints in the microenvironment of primary and metastatic renal cell cancer. Clin Cancer Res an Off J Am Assoc Cancer Res (2015) 21(13):3031–40. doi: 10.1158/1078-0432.Ccr-14-2926
7. McDermott D, Huseni M, Atkins M, Motzer R, Rini B, Escudier B, et al. Clinical activity and molecular correlates of response to atezolizumab alone or in combination with bevacizumab versus sunitinib in renal cell carcinoma. Nat Med (2018) 24(6):749–57. doi: 10.1038/s41591-018-0053-3
8. Helmink BA, Reddy SM, Gao J, Zhang S, Basar R, Thakur R, et al. B cells and tertiary lymphoid structures promote immunotherapy response. Nature (2020) 577(7791):549–55. doi: 10.1038/s41586-019-1922-8
9. Lozano AX, Chaudhuri AA, Nene A, Bacchiocchi A, Earland N, Vesely MD, et al. T Cell characteristics associated with toxicity to immune checkpoint blockade in patients with melanoma. Nat Med (2022) 28(2):353–62. doi: 10.1038/s41591-021-01623-z
10. Larkin J, Chiarion-Sileni V, Gonzalez R, Grob J, Cowey C, Lao C, et al. Combined nivolumab and ipilimumab or monotherapy in untreated melanoma. N Engl J Med (2015) 373(1):23–34. doi: 10.1056/NEJMoa1504030
11. Hellmann M, Callahan M, Awad M, Calvo E, Ascierto P, Atmaca A, et al. Tumor mutational burden and efficacy of nivolumab monotherapy and in combination with ipilimumab in small-cell lung cancer. Cancer Cell (2018) 33(5):853–61.e4. doi: 10.1016/j.ccell.2018.04.001
12. Au L, Hatipoglu E, Robert de Massy M, Litchfield K, Beattie G, Rowan A, et al. Determinants of anti-Pd-1 response and resistance in clear cell renal cell carcinoma. Cancer Cell (2021) 39(11):1497–518.e11. doi: 10.1016/j.ccell.2021.10.001
13. Braun DA, Street K, Burke KP, Cookmeyer DL, Denize T, Pedersen CB, et al. Progressive immune dysfunction with advancing disease stage in renal cell carcinoma. Cancer Cell (2021) 39(5):632–48.e8. doi: 10.1016/j.ccell.2021.02.013
14. Bi K, He MX, Bakouny Z, Kanodia A, Napolitano S, Wu J, et al. Tumor and immune reprogramming during immunotherapy in advanced renal cell carcinoma. Cancer Cell (2021) 39(5):649–61.e5. doi: 10.1016/j.ccell.2021.02.015
15. Gao J, Shi L, Zhao H, Chen J, Xiong L, He Q, et al. Loss of ifn-Γ pathway genes in tumor cells as a mechanism of resistance to anti-Ctla-4 therapy. Cell (2016) 167(2):397–404.e9. doi: 10.1016/j.cell.2016.08.069
16. Sade-Feldman M, Jiao Y, Chen J, Rooney M, Barzily-Rokni M, Eliane J, et al. Resistance to checkpoint blockade therapy through inactivation of antigen presentation. Nat Commun (2017) 8(1):1136. doi: 10.1038/s41467-017-01062-w
17. Shimizu S, Hiratsuka H, Koike K, Tsuchihashi K, Sonoda T, Ogi K, et al. Tumor-infiltrating Cd8 T-cell density is an independent prognostic marker for oral squamous cell carcinoma. Cancer Med (2019) 8(1):80–93. doi: 10.1002/cam4.1889
18. Litchfield K, Reading J, Lim E, Xu H, Liu P, Al-Bakir M, et al. Escape from nonsense-mediated decay associates with anti-tumor immunogenicity. Nat Commun (2020) 11(1):3800. doi: 10.1038/s41467-020-17526-5
19. Chowell D, Morris L, Grigg C, Weber J, Samstein R, Makarov V, et al. Patient hla class I genotype influences cancer response to checkpoint blockade immunotherapy. Sci (New York NY) (2018) 359(6375):582–7. doi: 10.1126/science.aao4572
20. Braun DA, Hou Y, Bakouny Z, Ficial M, Sant' Angelo M, Forman J, et al. Interplay of somatic alterations and immune infiltration modulates response to pd-1 blockade in advanced clear cell renal cell carcinoma. Nat Med (2020) 26(6):909–18. doi: 10.1038/s41591-020-0839-y
21. Eisenhauer E, Therasse P, Bogaerts J, Schwartz L, Sargent D, Ford R, et al. New response evaluation criteria in solid tumours: Revised recist guideline (Version 1.1). Eur J Cancer (2009) 45(2):228–47. doi: 10.1016/j.ejca.2008.10.026
22. Zou J, Deng F, Wang M, Zhang Z, Liu Z, Zhang X, et al. Sccode: An r package for data-specific differentially expressed gene detection on single-cell RNA-sequencing data. Brief Bioinf (2022). doi: 10.1093/bib/bbac180
23. Subramanian A, Tamayo P, Mootha V, Mukherjee S, Ebert B, Gillette M, et al. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA (2005) 102(43):15545–50. doi: 10.1073/pnas.0506580102
24. Jan B, Kosztyla D, von Törne C, Stenzinger A, Darb-Esfahani S, Dietel M, et al. Cancerclass: An r package for development and validation of diagnostic tests from high-dimensional molecular data. J Stat Software (2014) 59(1):1–19. doi: 10.18637/jss.v059.i01
25. Hao Y, Hao S, Andersen-Nissen E, Mauck WM 3rd, Zheng S, Butler A, et al. Integrated analysis of multimodal single-cell data. Cell (2021) 184(13):3573–87.e29. doi: 10.1016/j.cell.2021.04.048
26. Hänzelmann S, Castelo R, Guinney J. Gsva: Gene set variation analysis for microarray and RNA-seq data. BMC Bioinf (2013) 14:7. doi: 10.1186/1471-2105-14-7
27. Borcherding N, Vishwakarma A, Voigt AP, Bellizzi A, Kaplan J, Nepple K, et al. Mapping the immune environment in clear cell renal carcinoma by single-cell genomics. Commun Biol (2021) 4(1):122. doi: 10.1038/s42003-020-01625-6
28. Li H, van der Leun AM, Yofe I, Lubling Y, Gelbard-Solodkin D, van Akkooi ACJ, et al. Dysfunctional Cd8 T cells form a proliferative, dynamically regulated compartment within human melanoma. Cell (2019) 176(4):775–89.e18. doi: 10.1016/j.cell.2018.11.043
29. Truica M, Burns M, Han H, Abdulkadir S. Turning up the heat on myc: Progress in small-molecule inhibitors. Cancer Res (2021) 81(2):248–53. doi: 10.1158/0008-5472.Can-20-2959
30. Goel S, Decristo MJ, Watt AC, Brinjones H, Sceneay J, Li BB, et al. Cdk4/6 inhibition triggers anti-tumour immunity. Nature (2017) 548(7668):471–5. doi: 10.1038/nature23465
31. Kotanides H, Li Y, Malabunga M, Carpenito C, Eastman S, Shen Y, et al. Bispecific targeting of pd-1 and pd-L1 enhances T-cell activation and antitumor immunity. Cancer Immunol Res (2020) 8(10):1300–10. doi: 10.1158/2326-6066.Cir-20-0304
32. Katlinski K, Gui J, Katlinskaya Y, Ortiz A, Chakraborty R, Bhattacharya S, et al. Inactivation of interferon receptor promotes the establishment of immune privileged tumor microenvironment. Cancer Cell (2017) 31(2):194–207. doi: 10.1016/j.ccell.2017.01.004
33. Yang M, Li J, Gu P, Fan X. The application of nanoparticles in cancer immunotherapy: Targeting tumor microenvironment. Bioact Mater (2021) 6(7):1973–87. doi: 10.1016/j.bioactmat.2020.12.010
34. Cao Z, Kon N, Liu Y, Xu W, Wen J, Yao H, et al. An unexpected role for P53 in regulating cancer cell-intrinsic pd-1 by acetylation. Sci Adv (2021) 7(14):eabf4148. doi: 10.1126/sciadv.abf4148
35. Rahimi Kalateh Shah Mohammad G, Ghahremanloo A, Soltani A, Fathi E, Hashemy S. Cytokines as potential combination agents with pd-1/Pd-L1 blockade for cancer treatment. J Cell Physiol (2020) 235:5449–60. doi: 10.1002/jcp.29491
36. Deng L, Qian G, Zhang S, Zheng H, Fan S, Lesinski G, et al. Inhibition of mtor complex 1/P70 S6 kinase signaling elevates pd-L1 levels in human cancer cells through enhancing protein stabilization accompanied with enhanced Β-trcp degradation. Oncogene (2019) 38(35):6270–82. doi: 10.1038/s41388-019-0877-4
37. Wu J, Lu LY, Yu X. The role of Brca1 in DNA damage response. Protein Cell (2010) 1(2):117–23. doi: 10.1007/s13238-010-0010-5
38. Mahoney KM, Ross-Macdonald P, Yuan L, Song L, Veras E, Wind-Rotolo M, et al. Soluble pd-L1 as an early marker of progressive disease on nivolumab. J Immunother Cancer (2022) 10(2):e003527. doi: 10.1136/jitc-2021-003527
39. Jerby-Arnon L, Shah P, Cuoco MS, Rodman C, Su MJ, Melms JC, et al. A cancer cell program promotes T cell exclusion and resistance to checkpoint blockade. Cell (2018) 175(4):984–97.e24. doi: 10.1016/j.cell.2018.09.006
40. Grasso C, Tsoi J, Onyshchenko M, Abril-Rodriguez G, Ross-Macdonald P, Wind-Rotolo M, et al. Conserved interferon-Γ signaling drives clinical response to immune checkpoint blockade therapy in melanoma. Cancer Cell (2020) 38(4):500–15.e3. doi: 10.1016/j.ccell.2020.08.005
41. Johnson D, Estrada M, Salgado R, Sanchez V, Doxie D, Opalenik S, et al. Melanoma-specific mhc-ii expression represents a tumour-autonomous phenotype and predicts response to anti-Pd-1/Pd-L1 therapy. Nat Commun (2016) 7:10582. doi: 10.1038/ncomms10582
42. Roemer M, Redd R, Cader F, Pak C, Abdelrahman S, Ouyang J, et al. Major histocompatibility complex class ii and programmed death ligand 1 expression predict outcome after programmed death 1 blockade in classic Hodgkin lymphoma. J Clin Oncol Off J Am Soc Clin Oncol (2018) 36(10):942–50. doi: 10.1200/jco.2017.77.3994
43. Maugeri-Sacca M, Bartucci M, De Maria R. DNA Damage repair pathways in cancer stem cells. Mol Cancer Ther (2012) 11(8):1627–36. doi: 10.1158/1535-7163.MCT-11-1040
44. Motzer R, Rini B, McDermott D, Redman B, Kuzel T, Harrison M, et al. Nivolumab for metastatic renal cell carcinoma: Results of a randomized phase ii trial. J Clin Oncol Off J Am Soc Clin Oncol (2015) 33(13):1430–7. doi: 10.1200/jco.2014.59.0703
45. Choueiri T, Fishman M, Escudier B, McDermott D, Drake C, Kluger H, et al. Immunomodulatory activity of nivolumab in metastatic renal cell carcinoma. Clin Cancer Res an Off J Am Assoc Cancer Res (2016) 22(22):5461–71. doi: 10.1158/1078-0432.Ccr-15-2839
46. Miao D, Margolis C, Gao W, Voss M, Li W, Martini D, et al. Genomic correlates of response to immune checkpoint therapies in clear cell renal cell carcinoma. Sci (New York NY) (2018) 359(6377):801–6. doi: 10.1126/science.aan5951
47. Ascierto ML, McMiller TL, Berger AE, Danilova L, Anders RA, Netto GJ, et al. The intratumoral balance between metabolic and immunologic gene expression is associated with anti-Pd-1 response in patients with renal cell carcinoma. Cancer Immunol Res (2016) 4(9):726–33. doi: 10.1158/2326-6066.CIR-16-0072
48. Xiong D, Wang Y, You M. A gene expression signature of Trem2(Hi) macrophages and gammadelta T cells predicts immunotherapy response. Nat Commun (2020) 11(1):5084. doi: 10.1038/s41467-020-18546-x
49. Ayers M, Lunceford J, Nebozhyn M, Murphy E, Loboda A, Kaufman DR, et al. Ifn-Gamma-Related mrna profile predicts clinical response to pd-1 blockade. J Clin Invest (2017) 127(8):2930–40. doi: 10.1172/JCI91190
50. Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med (2018) 24(10):1550–8. doi: 10.1038/s41591-018-0136-1
51. Auslander N, Zhang G, Lee JS, Frederick DT, Miao B, Moll T, et al. Robust prediction of response to immune checkpoint blockade therapy in metastatic melanoma. Nat Med (2018) 24(10):1545–9. doi: 10.1038/s41591-018-0157-9
52. Yang S, Wu Y, Deng Y, Zhou L, Yang P, Zheng Y, et al. Identification of a prognostic immune signature for cervical cancer to predict survival and response to immune checkpoint inhibitors. Oncoimmunology (2019) 8(12):e1659094. doi: 10.1080/2162402X.2019.1659094
53. Dominguez CX, Muller S, Keerthivasan S, Koeppen H, Hung J, Gierke S, et al. Single-cell RNA sequencing reveals stromal evolution into Lrrc15(+) myofibroblasts as a determinant of patient response to cancer immunotherapy. Cancer Discov (2020) 10(2):232–53. doi: 10.1158/2159-8290.CD-19-0644
54. Huang AC, Orlowski RJ, Xu X, Mick R, George SM, Yan PK, et al. A single dose of neoadjuvant pd-1 blockade predicts clinical outcomes in resectable melanoma. Nat Med (2019) 25(3):454–61. doi: 10.1038/s41591-019-0357-y
55. Hugo W, Zaretsky JM, Sun L, Song C, Moreno BH, Hu-Lieskovan S, et al. Genomic and transcriptomic features of response to anti-Pd-1 therapy in metastatic melanoma. Cell (2016) 165(1):35–44. doi: 10.1016/j.cell.2016.02.065
56. Thompson JC, Hwang WT, Davis C, Deshpande C, Jeffries S, Rajpurohit Y, et al. Gene signatures of tumor inflammation and epithelial-to-Mesenchymal transition (Emt) predict responses to immune checkpoint blockade in lung cancer with high accuracy. Lung Cancer (2020) 139:1–8. doi: 10.1016/j.lungcan.2019.10.012
57. Friedlander P, Wassmann K, Christenfeld AM, Fisher D, Kyi C, Kirkwood JM, et al. Whole-blood RNA transcript-based models can predict clinical response in two Large independent clinical studies of patients with advanced melanoma treated with the checkpoint inhibitor, tremelimumab. J Immunother Cancer (2017) 5(1):67. doi: 10.1186/s40425-017-0272-z
58. Zhang Z, Wang Z, Chen Y, Wu H, Yin L, Zhao Q, et al. Integrated analysis of single-cell and bulk RNA sequencing data reveals a pan-cancer stemness signature predicting immunotherapy response. Genome Med (2022) 14(1):45. doi: 10.1186/s13073-022-01050-w
59. Braun D, Ishii Y, Walsh A, Van Allen E, Wu C, Shukla S, et al. Clinical validation of Pbrm1 alterations as a marker of immune checkpoint inhibitor response in renal cell carcinoma. JAMA Oncol (2019) 5(11):1631–3. doi: 10.1001/jamaoncol.2019.3158
60. Snyder A, Makarov V, Merghoub T, Yuan J, Zaretsky J, Desrichard A, et al. Genetic basis for clinical response to ctla-4 blockade in melanoma. New Engl J Med (2014) 371(23):2189–99. doi: 10.1056/NEJMoa1406498
61. Aggen D, Drake C. Biomarkers for immunotherapy in bladder cancer: A moving target. J Immunother Cancer (2017) 5(1):94. doi: 10.1186/s40425-017-0299-1
62. Fang W, Jin H, Zhou H, Hong S, Ma Y, Zhang Y, et al. Intratumoral heterogeneity as a predictive biomarker in anti-Pd-(L)1 therapies for non-small cell lung cancer. Mol Cancer (2021) 20(1):37. doi: 10.1186/s12943-021-01331-9
63. Hakimi AA, Voss MH, Kuo F, Sanchez A, Liu M, Nixon BG, et al. Transcriptomic profiling of the tumor microenvironment reveals distinct subgroups of clear cell renal cell cancer: Data from a randomized phase iii trial. Cancer Discov (2019) 9(4):510–25. doi: 10.1158/2159-8290.CD-18-0957
64. Choueiri TK, Albiges L, Haanen JBAG, Larkin JMG, Uemura M, Pal SK, et al. Biomarker analyses from javelin renal 101: Avelumab + axitinib (a+Ax) versus sunitinib (S) in advanced renal cell carcinoma (Arcc). J Clin Oncol (2019) 37(15_suppl):101. doi: 10.1200/JCO.2019.37.15_suppl.101
65. Lu S, Stein J, Rimm D, Wang D, Bell J, Johnson D, et al. Comparison of biomarker modalities for predicting response to pd-1/Pd-L1 checkpoint blockade: A systematic review and meta-analysis. JAMA Oncol (2019) 5(8):1195–204. doi: 10.1001/jamaoncol.2019.1549
66. Brahmer JR, Drake CG, Wollner I, Powderly JD, Picus J, Sharfman WH, et al. Phase I study of single-agent anti-programmed death-1 (Mdx-1106) in refractory solid tumors: Safety, clinical activity, pharmacodynamics, and immunologic correlates. J Clin Oncol (2010) 28(19):3167–75. doi: 10.1200/JCO.2009.26.7609
67. Dempke WCM, Fenchel K, Dale SP. Programmed cell death ligand-1 (Pd-L1) as a biomarker for non-small cell lung cancer (Nsclc) treatment-are we barking up the wrong tree? Transl Lung Cancer Res (2018) 7(Suppl 3):S275–S9. doi: 10.21037/tlcr.2018.04.18
68. Borst J, Ahrends T, Babala N, Melief CJM, Kastenmuller W. Cd4(+) T cell help in cancer immunology and immunotherapy. Nat Rev Immunol (2018) 18(10):635–47. doi: 10.1038/s41577-018-0044-0
69. Melssen M, Slingluff C. Vaccines targeting helper T cells for cancer immunotherapy. Curr Opin Immunol (2017) 47:85–92. doi: 10.1016/j.coi.2017.07.004
70. Castellino F, Germain R. Cooperation between Cd4+ and Cd8+ T cells: When, where, and how. Annu Rev Immunol (2006) 24:519–40. doi: 10.1146/annurev.immunol.23.021704.115825
71. Savage P, Klawon D, Miller C. Regulatory T cell development. Annu Rev Immunol (2020) 38:421–53. doi: 10.1146/annurev-immunol-100219-020937
72. Tay RE, Richardson EK, Toh HC. Revisiting the role of Cd4(+) T cells in cancer immunotherapy-new insights into old paradigms. Cancer Gene Ther (2021) 28(1-2):5–17. doi: 10.1038/s41417-020-0183-x
73. Kamada T, Togashi Y, Tay C, Ha D, Sasaki A, Nakamura Y, et al. Pd-1(+) regulatory T cells amplified by pd-1 blockade promote hyperprogression of cancer. Proc Natl Acad Sci USA (2019) 116(20):9999–10008. doi: 10.1073/pnas.1822001116
74. Matoba T, Imai M, Ohkura N, Kawakita D, Ijichi K, Toyama T, et al. Regulatory T cells expressing abundant ctla-4 on the cell surface with a proliferative gene profile are key features of human head and neck cancer. Int J Cancer (2019) 144(11):2811–22. doi: 10.1002/ijc.32024
75. Sanmamed MF, Nie X, Desai SS, Villaroel-Espindola F, Badri T, Zhao D, et al. A burned-out Cd8(+) T-cell subset expands in the tumor microenvironment and curbs cancer immunotherapy. Cancer Discov (2021) 11(7):1700–15. doi: 10.1158/2159-8290.CD-20-0962
76. Kamphorst AO, Pillai RN, Yang S, Nasti TH, Akondy RS, Wieland A, et al. Proliferation of pd-1+ Cd8 T cells in peripheral blood after pd-1-Targeted therapy in lung cancer patients. Proc Natl Acad Sci USA (2017) 114(19):4993–8. doi: 10.1073/pnas.1705327114
77. Kurtulus S, Madi A, Escobar G, Klapholz M, Nyman J, Christian E, et al. Checkpoint blockade immunotherapy induces dynamic changes in pd-1(-)Cd8(+) tumor-infiltrating T cells. Immunity (2019) 50(1):181–94.e6. doi: 10.1016/j.immuni.2018.11.014
78. Dai S, Zeng H, Liu Z, Jin K, Jiang W, Wang Z, et al. Intratumoral Cxcl13(+)Cd8(+)T cell infiltration determines poor clinical outcomes and immunoevasive contexture in patients with clear cell renal cell carcinoma. J Immunother Cancer (2021) 9(2):e001823. doi: 10.1136/jitc-2020-001823
79. Zheng Y, Chen Z, Han Y, Han L, Zou X, Zhou B, et al. Immune suppressive landscape in the human esophageal squamous cell carcinoma microenvironment. Nat Commun (2020) 11(1):6268. doi: 10.1038/s41467-020-20019-0
80. Im SJ, Hashimoto M, Gerner MY, Lee J, Kissick HT, Burger MC, et al. Defining Cd8+ T cells that provide the proliferative burst after pd-1 therapy. Nature (2016) 537(7620):417–21. doi: 10.1038/nature19330
81. Braun NA, Celada LJ, Herazo-Maya JD, Abraham S, Shaginurova G, Sevin CM, et al. Blockade of the programmed death-1 pathway restores sarcoidosis Cd4(+) T-cell proliferative capacity. Am J Respir Crit Care Med (2014) 190(5):560–71. doi: 10.1164/rccm.201401-0188OC
82. Ugel S, Cane S, De Sanctis F, Bronte V. Monocytes in the tumor microenvironment. Annu Rev Pathol (2021) 16:93–122. doi: 10.1146/annurev-pathmechdis-012418-013058
83. Jakubzick C, Randolph G, Henson P. Monocyte differentiation and antigen-presenting functions. Nat Rev Immunol (2017) 17(6):349–62. doi: 10.1038/nri.2017.28
84. González-Domínguez É, Samaniego R, Flores-Sevilla J, Campos-Campos S, Gómez-Campos G, Salas A, et al. Cd163l1 and Clec5a discriminate subsets of human resident and inflammatory macrophages in vivo. J leukocyte Biol (2015) 98(4):453–66. doi: 10.1189/jlb.3HI1114-531R
85. Tang-Huau T, Gueguen P, Goudot C, Durand M, Bohec M, Baulande S, et al. Human in vivo-generated monocyte-derived dendritic cells and macrophages cross-present antigens through a vacuolar pathway. Nat Commun (2018) 9(1):2570. doi: 10.1038/s41467-018-04985-0
86. Lam KC, Araya RE, Huang A, Chen Q, Di Modica M, Rodrigues RR, et al. Microbiota triggers sting-type I ifn-dependent monocyte reprogramming of the tumor microenvironment. Cell (2021) 184(21):5338–56.e21. doi: 10.1016/j.cell.2021.09.019
87. Haghverdi L, Lun ATL, Morgan MD, Marioni JC. Batch effects in single-cell RNA-sequencing data are corrected by matching mutual nearest neighbors. Nat Biotechnol (2018) 36(5):421–7. doi: 10.1038/nbt.4091
88. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res (2015) 43(7):e47. doi: 10.1093/nar/gkv007
Keywords: clear cell renal cell carcinoma (ccRCC), immune checkpoint blockade (ICB) therapy, response prediction model, single-cell RNA-seq, machine learning (ML)
Citation: Zheng K, Gao L, Hao J, Zou X and Hu X (2022) An immunotherapy response prediction model derived from proliferative CD4+ T cells and antigen-presenting monocytes in ccRCC. Front. Immunol. 13:972227. doi: 10.3389/fimmu.2022.972227
Received: 18 June 2022; Accepted: 08 August 2022;
Published: 25 August 2022.
Edited by:
Lewis Zhichang Shi, University of Alabama at Birmingham, United StatesReviewed by:
Ciputra Hartana, Ragon Institute, United StatesYuming Yu, Guangdong Provincial People’s Hospital, China
Copyright © 2022 Zheng, Gao, Hao, Zou and Hu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Jie Hao, jhao@fudan.edu.cn; Xin Zou, xzou@fudan.edu.cn; Xiaoyong Hu, huxiaoyong@sjtu.edu.cn
†These authors have contributed equally to this work