Integrated analysis of single-cell and bulk RNA sequencing data reveals a pan-cancer stemness signature predicting immunotherapy response
Genome Medicine volume 14, Article number: 45 (2022)
Although immune checkpoint inhibitor (ICI) is regarded as a breakthrough in cancer therapy, only a limited fraction of patients benefit from it. Cancer stemness can be the potential culprit in ICI resistance, but direct clinical evidence is lacking.
Publicly available scRNA-Seq datasets derived from ICI-treated patients were collected and analyzed to elucidate the association between cancer stemness and ICI response. A novel stemness signature (Stem.Sig) was developed and validated using large-scale pan-cancer data, including 34 scRNA-Seq datasets, The Cancer Genome Atlas (TCGA) pan-cancer cohort, and 10 ICI transcriptomic cohorts. The therapeutic value of Stem.Sig genes was further explored using 17 CRISPR datasets that screened potential immunotherapy targets.
Cancer stemness, as evaluated by CytoTRACE, was found to be significantly associated with ICI resistance in melanoma and basal cell carcinoma (both P < 0.001). Significantly negative association was found between Stem.Sig and anti-tumor immunity, while positive correlations were detected between Stem.Sig and intra-tumoral heterogenicity (ITH) / total mutational burden (TMB). Based on this signature, machine learning model predicted ICI response with an AUC of 0.71 in both validation and testing set. Remarkably, compared with previous well-established signatures, Stem.Sig achieved better predictive performance across multiple cancers. Moreover, we generated a gene list ranked by the average effect of each gene to enhance tumor immune response after genetic knockout across different CRISPR datasets. Then we matched Stem.Sig to this gene list and found Stem.Sig significantly enriched 3% top-ranked genes from the list (P = 0.03), including EMC3, BECN1, VPS35, PCBP2, VPS29, PSMF1, GCLC, KXD1, SPRR1B, PTMA, YBX1, CYP27B1, NACA, PPP1CA, TCEB2, PIGC, NR0B2, PEX13, SERF2, and ZBTB43, which were potential therapeutic targets.
We revealed a robust link between cancer stemness and immunotherapy resistance and developed a promising signature, Stem.Sig, which showed increased performance in comparison to other signatures regarding ICI response prediction. This signature could serve as a competitive tool for patient selection of immunotherapy. Meanwhile, our study potentially paves the way for overcoming immune resistance by targeting stemness-associated genes.
Immune checkpoint inhibitor (ICI) has ushered in a new era of cancer treatment and provided unprecedented clinical benefits for patients . However, only a relatively small proportion of patients respond to it , which highlights the necessity of biomarker research for optimizing patient selection and combination strategies to tackle immune resistance.
Traditional biomarker research mostly focused on the analysis of whole exome sequencing (WES) or RNA sequencing (RNA-Seq) from intact tumor tissue (bulk data) [3,4,5,6,7,8], which only reflects the average genetic profile across a large population of different cells. Pre-existing ICI biomarkers derived from these studies showed limited predictive values. The development of single-cell RNA sequencing (scRNA-Seq) enables us to dissect gene expression at single-cell resolution and identify novel biomarkers with better performance .
Cancer stem cells (CSCs) are self-renewal cells that promote tumor initiation, progression, and metastasis . Mounting evidences revealed a prominent association between stemness and cancer immune evasion and resistance . A previous study demonstrated that high stemness correlates with immune cell exclusion across 21 solid cancer types , but direct clinical evidence validating the negative association between stemness and ICI outcomes is lacking. With the help of a powerful computational framework (CytoTRACE) developed by Gulati et al., we can accurately characterize cancer stemness and identify stemness-correlated genes at the resolution of single-cell level to better investigate the impacts of stemness on ICI .
In this study, we revealed and verified the negative association between cancer stemness and ICI outcomes in two scRNA-Seq ICI cohorts [14, 15]. Thereafter a stemness signature (Stem.Sig) was developed through an integrative analysis of 34 scRNA-Seq datasets, which consisted of 345 patients and 663,760 cells across 17 cancer types [14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43]. The predictive value of Stem.Sig was further explored and validated through a comprehensive analysis of pan-cancer transcriptomic data (10,154 patients; 30 cancer types) , 17 CRISPR datasets (4 cancer types) [45,46,47,48,49,50,51], and 10 independent ICI cohorts (921 patients; 5 cancer types) [52,53,54,55,56,57,58,59,60,61]. Our findings uncovered the potential of Stem.Sig for predicting ICI outcomes more accurately than previously recognized signatures across multiple cancer types.
scRNA-Seq ICI cohorts
To investigate the relationship between cancer cell stemness and immunotherapy efficacy, a melanoma cohort with both ICI response and scRNA-Seq data was analyzed [14, 15]. Another independent scRNA-seq ICI cohort of basal cell carcinoma was used to validate the results. Data of these two cohorts was accessed through GEO accession number: GSE115978  and GSE123813 , respectively (Additional file 1: Table S1).
Pan-cancer scRNA-Seq datasets
For the development of stemness signature (Stem.Sig), 34 scRNA-Seq datasets with both malignant and stromal/immune cells data were collected from the TISCH portal (http://tisch.comp-genomics.org/) , which consist of 345 patients and 663,760 cells (Additional file 1: Table S2). These datasets cover 17 cancer types, including basal cell carcinoma (BCC), breast cancer (BRCA), cholangiocarcinoma (CHOL), colorectal cancer (CRC), glioma, head and neck cancer (HNSC), liver hepatocellular carcinoma (LIHC), medulloblastoma (MB), Merkel cell carcinoma (MCC), multiple myeloma (MM), neuroendocrine tumor (NET), non-small cell lung cancer (NSCLC), ovarian serous cystadenocarcinoma (OV), pancreatic adenocarcinoma (PAAD), skin cutaneous melanoma (SKCM), stomach adenocarcinoma (STAD), and uveal melanoma (UVM) [14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43].
Pan-cancer transcriptomic data
Transcriptomic data of The Cancer Genome Atlas (TCGA) Pan-cancer cohort was downloaded from the UCSC Xena data portal (https://xenabrowser.net)  to explore the potential links between Stem.Sig and immune suppression across 30 different cancer types. Three cancer types were excluded from our analysis, including diffuse large B cell lymphoma (DLBC), acute myeloid leukemia (LAML), and thymoma (THYM), as they mainly consist of immune cells . Total mutation burden (TMB) was retrieved from cBioPortal (https://www.cbioportal.org) [64, 65] and intratumor heterogeneity (ITH) data was from Thorsson et al. , which were used for analyzing the correlation between Stem.Sig and TMB or ITH.
ICI RNA-Seq cohorts
To validate the predictive value of Stem.Sig, we systemically collected transcriptomic data and clinical information of pretreatment samples from 10 ICI RNA-Seq cohorts, including 5 SKCM cohorts (Hugo 2016 , Liu 2019 , Gide 2019 , Riaz 2017 , Van Allen 2015 ), 2 urothelial carcinoma (UC) cohorts (Mariathasan 2018 , Synder 2017 ), 1 glioblastoma multiform (GBM) cohort (Zhao 2019 ), 1 gastric cancer (GC) cohort (Kim 2018 ) and 1 renal cell carcinoma (RCC) cohort (Braun 2020 ). Anti-PD-1 therapy, anti-PD-L1 therapy, anti-CTLA4 therapy, and anti-PD-(L)1 plus anti-CTLA-4 combination therapy were employed in 6, 2, 1, and 1 cohort, respectively. Cohort Hugo 2016  comprises 27 pre-treated tumor samples from 26 patients, while cohort Zhao 2019  consists of 34 pre-treated tumor samples from 17 patients. In these two cohorts, we randomly selected a single tumor sample for each corresponding patient. The details of these cohorts are summarized in Additional file 1: Table S3.
CRISPR screening data
To explore potential therapeutic targets of Stem.Sig genes, we collected data from 7 published CRISPR/Cas9 screening studies that assessed the individual effect of each gene knockout on tumor immunity, including Freeman 2019 , Kearney 2018 , Manguso 2017 , Pan 2018 , Patel 2017 , Vredevoogd 2019 , and Lawson 2020 . The first six CRISPR studies have been previously curated by Fu et al.  In addition to Fu et al., we further collected another CRISPR cohort from Lawson et al. . According to model cell lines and treatment conditions applied, these seven CRISPR studies were divided into 17 datasets (Additional file 1: Table S4). The CRISPR analysis covers melanoma, breast cancer, colon cancer, and renal cancer cell lines. We utilized these data to identify genes that are more likely to modulate lymphocyte-mediated cancer killing and influence immunotherapy response across different datasets.
The process of CRISPR screens is to perform genome-wide CRISPR-Cas9 knockout across various cancer cell lines that were co-cultured with/without cytotoxic lymphocytes (CTLs) in vitro or implanted into immune-deficient mice/immune-competent mice in vivo. Then RNA sequencing is used to estimate the abundance of sgRNA targeting the corresponding gene. To measure the effect of gene knockout on cancer fitness under the pressure of CTLs or anti-tumor immunity, log-fold changes of sgRNA reads are calculated for paired screens of cell lines (with CTLs vs. without CTLs; immune-deficient mice vs. immune-competent mice) . Normalized z scores were called from the log-fold changes in order to remove batch effects and compare genes among CRISPR datasets from different studies. The lower z scores indicate better immune response after gene knockout. We also ranked genes based on the average z scores across 17 datasets. Top-ranked genes with lower z scores were characterized as immune resistant.
scRNA-Seq data analysis
Previous studies revealed a global reduction of chromatin accessibility during lineage commitment . Since chromatin accessibility could be quantitatively reflected by single-cell gene counts, Gulati et al. had discovered a prominent association between single-cell gene counts and differential status of the corresponding cell . Higher single-cell gene counts correlate with less cellular differentiation (higher stemness). The CytoTRACE algorithm is developed by Gulati et al. to capture, smooth, and calculate the expression level of genes that are most highly correlated with single-cell gene counts with scRNA-Seq data. When the calculation of the CytoTRACE algorithm is finished, each single cell will get a score that represents its stemness within the given dataset. CytoTRACE is a robust computational framework for differentiation states prediction via scRNA-seq data, which was validated in large-scale datasets and outperformed pre-existing computational techniques of stemness . R package CytoTRACE v0.3.3 was applied to calculate the CytoTRACE scores for malignant cells. CytoTRACE scores range from 0 to 1, while higher scores indicate higher stemness (less differentiation) and vice versa.
The R package Seurat v4.0.6 was used to identify differentially expressed genes of malignant cells in each dataset. Genes with the log-fold change (logFC) ≥ 0.25 and false discovery rate (FDR) < 1e−05 were considered as differentially up-regulated genes of malignant cells in each dataset .
Anti-tumor immunity and pathway analysis
Over Representation Analysis (ORA) was conducted to determine whether known biological functions or processes are enriched in Stem.Sig . The R package clusterProfiler v4.2.1 was used to perform ORA and visualize the results .
Gene set variation analysis (GSVA) was used to calculate the Stem.Sig scores and relative enrichment of other gene signatures and biological pathways across sample populations. The R package GSVA v1.42.0 was applied to perform GSVA .
We further evaluated the correlation between Stem.Sig and tumor-infiltrating leukocytes (TILs)/immune-related genes (IRGs) in TCGA cohort. Immune-related genes and their functional classifications were obtained from Thorsson et al. . The R package MCP-counter v1.1.0 was utilized to estimate the abundance of tumor-infiltrating leukocytes .
The primary clinical outcomes were objective response rate (ORR) and overall survival (OS). ORR was assessed using Response Evaluation Criteria in Solid Tumors (RECIST) version 1.1 in all cohorts , except cohort Hugo 2016 , whose ORR was assessed using immune-related RECIST (irRECIST). Patients were divided into two groups according to their response status: complete response (CR) and partial response (PR) as responders, or stable disease (SD) and progressive disease (PD) as non-responders.
Derivation of predictive model for ICI response
Top five ICI RNA-Seq cohorts with most patients were combined to form a large cohort (n= 772), including RCC (n=181), UC (n=348), and SKCM (n=243). These five cohorts are Braun 2020 RCC , Mariathasan 2018 UC , Liu 2019 SKCM , Gide 2019 SKCM , and Riaz 2017 SKCM . We used ComBat method to remove the batch effect of different ICI RNA-Seq cohorts . Then we randomly split this combined cohort into two datasets: training set (80%, n = 618) and validation set (20%, n=154). The other five ICI RNA-Seq cohorts [57,58,59,60,61] were consolidated as an independent testing set (n = 149).
Model training and parameter tuning
We trained the ICI response classification model with Stem.Sig, using seven common machine learning (ML) algorithms, including support vector machine (SVM), Naïve Bayes (“NB”), random forest (“RF”), k-nearest neighbors (“KNN”), AdaBoost Classification Trees (“AdaBoost”), boosted logistic regressions (“LogiBoost”), and cancerclass [73, 74]. For each ML algorithm with parameters except cancerclass, fivefold cross-validation (CV) was adopted for hyperparameter tuning to optimize the performance of the model. To ensure robustness, we repeated the optimization process 10 times with different random seeds for each single resampling . As for cancerclass which does not require parameters, we trained the model using the entire training set directly.
Model validation and independent testing
We had seven models derived from the training set, using different ML algorithms. Then we applied these models to the validation set and compared their results. The model with the best performance was chosen as the final Stem.Sig model. To evaluate the predictive value of the final model, we applied it to the testing set.
Comparing Stem.Sig with other predictive gene signatures
To further evaluate the predictive value of Stem.Sig, we compared Stem.Sig with other ICI response signatures reported previously, including six pan-cancer signatures (INFG.Sig , T.cell.inflamed.Sig , PDL1.Sig , LRRC15.CAF.Sig , NLRP3.Sig , and Cytotoxic.Sig ) and seven melanoma-specific signatures (CRMA.Sig , IMPRES.Sig , IPRES.Sig , TcellExc.Sig , ImmmunCells.Sig , IMS.Sig , and TRS.Sig ). Pan-cancer signatures were compared with Stem.Sig in the testing set regarding the performance of ICI response prediction. As for melanoma-specific signatures, we compared their performance with Stem.Sig using melanoma patients from the testing set (Hugo 2016  and Van Allen 2015 ). Codes and algorithms for the 13 aforementioned signatures were derived from their original studies, such as ssGSEA for NLRP3.Sig , cancerclass for ImmuneCell.Sig , overall expression for TcellExc.Sig , and so on. Details of these signatures and their corresponding algorithms can be found in Additional file 1: Table S5.
Statistical analyses were performed using R v4.1.1 (https://www.r-project.org). Comparison of CytoTRACE scores between response and non-response subgroups was analyzed by two-sided Wilcoxon tests. We used Spearman correlation to evaluate the association between Stem.Sig and biological pathways or immune features. Benjamini-Hochberg procedure (B-H) was applied to calculate FDR. The R package caret v6.0-90 and cancerclass v1.34.0 were used for model training, validation, and testing . The receiver operating characteristic (ROC) curve was used and a larger area under the ROC curve (AUC) indicated a better predictive performance An AUC of 0.9–1.0 is considered excellent, 0.8–0.9 very good, 0.7–0.8 good, 0.6–0.7 sufficient, 0.5–0.6 bad, and less than 0.5 considered not useful . Patients predicted by the final model as “NR” and “R” were categorized into high-risk and low-risk subgroups for survival analysis. The association between the Stem.Sig-based risk grouping and OS was analyzed by Cox proportional hazards regression analysis. Further survival analysis of individual cohort from testing set was adjusted for available confounding factors, including TMB, tumor purity, sex, and age.
Cancer stemness is associated with ICI resistance
A previous published ICI SKCM cohort with scRNA-seq data was firstly employed to evaluate the association between cancer stemness and ICI outcomes . After removing patients without malignant cells data, we adopted 24 patients from this cohort, consisting of 11 non-responders (NR) and 13 treatment-naïve patients (TN). Ideally, it is better to compare the cancer stemness between responders (R) and non-responders. However, data of responders was not available in this cohort. Given that treatment naïve patients likely include both potential responders and non-responders, comparison of stemness was conducted between NR and TN as previously described . As shown in Fig. 1A, cancer cells with high stemness were enriched in the NR subgroup. Further analysis showed that tumors from the NR subgroup had a significantly higher level of stemness (P < 0.001, Fig. 1B), indicating that cancer stemness is negatively associated with ICI outcomes. Another ICI cohort with a different cancer type (BCC) was employed to validate this finding . In the BCC cohort, tumor stemness of 4 non-responders was compared to that of 6 responders. We found a more prominent gap of stemness level between NR and R subgroups in the BCC cohort (P < 0.001, Fig. 1C and D).
Development of Stem.Sig through pan-cancer scRNA analysis
As cancer stemness is significantly associated with ICI resistance, we hypothesized that a Stem.Sig reflecting the stemness level of the tumor may help in the prediction of ICI efficacy. Therefore, 34 scRNA-Seq datasets were employed to develop the Stem.Sig (Fig. 2A; Additional file 1: Table S6). We performed Spearman correlation analysis between gene expression level and CytoTRACE scores for malignant cells among pan-cancer scRNA datasets. Genes that were positively correlated with CytoTRACE scores (Spearman R > 0 and FDR < 1e−05) were regarded as Gx. Genes that were differentially up-regulated in malignant cells were regarded as Gy. To obtain up-regulated tumor-specific genes that were positively associated with stemness, Gx and Gy were intersected to give rise to Gn for each dataset . For example, G1 consisted of genes derived from the intersection of Gx and Gy in the first scRNA-Seq dataset. Geometric mean of Spearman R was calculated for each gene across G1–G34. Finally, genes with geometric mean of Spearman R > 0.4 (moderate to strong correlation) were pooled as Stem.Sig .
We investigated the biological functions that were over-represented in Stem.Sig (Fig. 2B). The enriched pathways mainly comprise processes involving hypoxia, glycolysis, ubiquitination, EPH-ephrin signaling, WNT signaling, and nucleotide excision repair (NER). Specific genes of these pathways were shown in the cnetplot of Fig. 2B. Some genes have been reported to be associated with unfavorable outcomes of immunotherapy, such as EPHA3, EPHA7, ENO1, ACTG1, DKK2, NPM1, and BCL10 [6, 88,89,90,91,92].
Analysis of the potential links between Stem.Sig and immune suppression using pan-cancer TCGA cohort
First, we performed a thorough analysis of Stem.Sig and 75 immune-related genes . A general negative association was observed between Stem.Sig and expression level of immune-related genes across 30 different cancer types (Fig. 3A). Then we evaluated the infiltration status of immune cells to better characterize the tumor immune microenvironment (TIME). Tumors with high Stem.Sig had decreased cytotoxic immune cells, including CD8+ T cells, NK cells and macrophages (Fig. 3B). Taken together, these results indicated that Stem.Sig was negatively associated with anti-tumor immunity.
Secondly, we analyzed the enrichment of hallmark pathways regarding the expression level of Stem.Sig to investigate whether immunosuppressive biological functions were upregulated in high Stem.Sig tumors. Metabolic pathways, DNA repair, and MYC signaling were found to be enriched in tumors with high Stem.Sig(Fig. 3C). All of these pathways contributed to the poor immune response according to previous studies [93,94,95].
Furthermore, we investigated the association between Stem.Sig and ITH, a stemness-associated feature that mediates immunosuppression . As expected, ITH were positively correlated with Stem.Sig across cancer types. (R = 0.42, P = 0.021, Fig. 3D). Same analysis was applied to TMB, a well-known immune-relevant factor. Interestingly, a similar positive association was detected between Stem.Sig and TMB as well (R = 0.47, P = 0.008, Fig. 3E). High TMB indicates better immune response, while high Stem.Sig dose the opposite. To better elucidate the correlation of anti-tumor immunity with both Stem.Sig and TMB, we further divide patients into four subgroups: high Stem.Sig/high TMB (HSHT), high Stem.Sig/low TMB (HSLT), low Stem.Sig / high TMB (LSHT), and low Stem.Sig / low TMB (LSLT). Median GSVA score of Stem.Sig and median TMB were used as thresholds for grouping. Then we compared the abundance of immune cells among these four subgroups. Interestingly, LSHT was found with the highest level of cytotoxic lymphocytes (P < 0.001), while HSLT with the lowest level (P<0.001). High cancer stemness (HS) could promote immune resistance and evasion, while low TMB level (LT) is associated with decreased anti-tumor immunity due to lack of antigenicity. As expected, decreased infiltration of cytotoxic lymphocytes was detected for both HS and LT groups (p < 0.001, Additional file 2: Fig.S1 A and B). It is reasonable that the coexistence of these two factors (HSLT) may result in a TIME with the least infiltration of cytotoxic lymphocytes. On the contrary, LSHT could lead to the most abundant CTLs in the TIME. However, the anti-tumor immunity of the other two groups (HSHT and LSLT) seems to be more controversial than the aforementioned groups (HSLT, LSHT), since HSHT and LSLT both have an immune-suppressed (HS or LT) factor and an immune-promoted (LS or HT) factor. Further subgroup analysis found a higher level of cytotoxic lymphocytes in LSLT than in HSHT (p < 0.001, Additional file 2: Fig.S1 C). In conclusion, the order of anti-immunity from highest to lowest is: LSHT > LSLT > HSHT > HSLT (all p < 0.001, Additional file 2: Fig.S1 C). Therefore, tumors with low Stem.Sig presented with significantly better anti-tumor immunity than those with high Stem.Sig regardless of TMB level.
Immunotherapy outcome prediction by Stem.Sig
To investigate the predictive value of Stem.Sig, we collected bulk RNA-Seq data and clinical information from 10 ICI cohorts. Pre-treatment samples of these cohorts were curated and analyzed. Patients received anti-PD(L)-1, anti-CTLA-4, or anti-PD(L)-1 plus anti-CTLA-4. All these 10 cohorts were split into 3 data set: training set (n=620), validation set (n=154), and testing set (n=149). The flow chart of the analysis process was shown in Fig. 4A. Firstly, we trained the model with seven different machine learning algorithms and applied 10-time repeated 5-fold cross-validation for parameter optimization of each model. After training, we harvested seven models. Then, we evaluated and compared the AUC of these models in the validation cohort. Naïve Bayes model achieved the highest AUC of 0.71 and was selected as Stem.Sig model (Fig. 4B). For further assessment of the Stem.Sig model, we applied it to the independent testing set to predict ICI response and observed a same AUC of 0.71 (Fig.4C).
To evaluate whether the Stem.Sig model can predict overall survival, we divided ICI-treated patients into low-risk and high-risk subgroups based on the predicted “R” and “NR” respectively. The Kaplan-Meier analysis of OS was shown in Fig. 4D. Low-risk group achieved a significantly longer overall survival in training, validation, and testing sets (all log-rank p < 0.01). In the validation cohort, high-risk patients predicted by the Stem.Sig model had a median OS of only 13.3 months, compared to 31.2 months of low-risk patients (HR: 1.87; 95%CI: 1.21–2.90). In the testing set, a similar median OS of 13.4 months was observed in high-risk patients, while low-risk ones had not reached the median OS (HR: 3.08; 95%CI: 1.64–5.81).
We performed subgroup analysis for five individual cohorts that contribute to the testing set. Regarding ICI response prediction, AUC ranged from 0.62 to 0.81 among these cohorts (Additional file 2: Fig.S2A). Van Allen 2015 SKCM achieved a favorable AUC of 0.81 (95%CI: 0.66−0.95), followed by Synder 2017 UC (AUC: 0.80; 95%CI: 0.61−0.99). Compared to other cohorts, Zhao 2019 GBM presented with the lowest AUC of 0.62 (95%CI: 0.33−0.91). In survival analysis, Kim 2018 GC was removed due to a lack of OS data. For the other four cohorts, we observed a HR ranged from 1.73 to 4.05 in high-risk patients predicted by the Stem.Sig model (Additional file 2: Fig.S2B). After adjusting available confounding factors, significant survival benefits were still found in Van Allen 2015 SKCM (adjusted p = 0.02) and Synder 2017 UC (adjusted p = 0.02), while the other two cohorts showed only numerical survival differences. It is possibly due to the limited sample size.
We further compared the performance of Stem.Sig with previous well-established predictive gene signatures. Compared with pan-cancer signatures (INFG.Sig , T.cell.inflamed.Sig , PDL1.Sig , LRRC15.CAF.Sig , NLRP3.Sig , and Cytotoxic.Sig ), Stem.Sig showed best performance in the testing set with an AUC of 0.71, followed by INFG.Sig with an AUC of 0.66 (Fig. 5A). Most of these pan-cancer signatures showed ideal performance in only one or two cohorts. For example, AUC of INFG.Sig reached 0.85 in Kim 2018 GC and 0.67 in Van Allen 2015 SKCM, but it decreased to 0.53–0.54 in the other three cohorts (Additional file 1: Table S7). However, Stem.Sig achieved sufficient to very good performance in all cohorts, covering four cancer types: SKCM, GBM, UC, and GC, which further stresses its potential as a predictive model of ICI response in a pan-cancer manner (Fig. 5B). Compared with melanoma-specific signatures (CRMA.Sig , IMPRES.Sig , IPRES.Sig , TcellExc.Sig , ImmmunCells.Sig , IMS.Sig , and TRS.Sig ), Stem.Sig remained in the top 3 with an AUC of 0.76 in prediction of ICI response regarding melanoma patients. IMPRES.Sig and CRMA.Sig showed a slightly better AUC of 0.81 and 0.77 than Stem.Sig.
Exploration of potential therapeutic targets from Stem.Sig using CRISPR screen data
We systemically collected immune response data of knockout genes from seven CRISPR cohorts, which were further divided into 17 datasets according to the model cells and treatment conditions used in these CRISPR cohorts. Totally, there were 22,505 genes recorded by these CRISPR datasets. We ranked genes based on their mean z scores. Top-ranked genes were immune-resistant genes, which may promote anti-tumor immunity after knockout. Bottom-ranked genes were immune-sensitive genes, which may suppress anti-tumor immunity after knockout. The process of gene ranking was shown in Fig. 6A. Among all 22,505 genes, the number of 1%, 2%, and 3% top-ranked genes was 225, 450, and 675, respectively. Next, we calculated the percentage of top-ranked genes that were presented in Stem.Sig and previous immune-resistant signatures, including TcellExc.Sig, ImmuneCells.Sig, IMS.Sig, LRRC15.CAF.Sig, and CRMA.Sig (except IPRES.Sig, which comprises 73 genetic pathways instead of individual genes) [14, 78, 81, 83, 84]. Stem.Sig, TcellExc.Sig, IMS.Sig, and ImmuneCells.Sig were the only four gene sets that had genes ranked in the top 3%. As expected, Stem.Sig had the highest percentage of top-ranked genes than other signatures (Fig. 6B). Immune-resistant genes (3% top-ranked genes) were significantly over-represented in Stem.Sig (P=0.03; Fisher’s exact test). There were 20 genes of Stem.Sig that were ranked in the top 3%, including EMC3, BECN1, VPS35, PCBP2, VPS29, PSMF1, GCLC, KXD1, SPRR1B, PTMA, YBX1, CYP27B1, NACA, PPP1CA, TCEB2, PIGC, NR0B2, PEX13, SERF2, and ZBTB43. Immune-resistant features of these stemness-associated genes were validated by multiple independent CRISPR datasets (Fig. 6C), which may serve as potential therapeutic targets in synergy with ICB.
Although the mechanism between cancer stemness and anti-tumor immunity has been widely explored [10, 12, 96, 97], direct clinical evidence on the association of stemness and ICI response has not been reported. Here we utilized CytoTRACE to evaluate the stemness level of individual malignant cells and uncovered the inverse correlation between stemness and ICI outcomes, supported by the results from two ICI scRNA-Seq cohorts of SKCM and BCC [14, 15]. CSCs have been found in virtually all solid tumors . Motivated by these observations, we hypothesized that the negative association between stemness and ICI efficacy generally existed across various cancers. Therefore, a large-scale comprehensive analysis was performed to identify over-expressed genes in malignant cells that significantly correlated with increased stemness. These genes formed a pan-cancer stemness signature, namely, Stem.Sig. We carefully validated the predictive value of Stem.Sig. Remarkably, Stem.Sig achieved better performance of predicting ICI response than previous predictive signatures across multiple independent ICI cohorts with bulk RNA-Seq data [52,53,54,55,56,57,58,59,60,61]. This study is the first report to demonstrate the robust link between stemness and ICI outcomes through a comprehensive analysis of large-scale data. Most importantly, we constructed a gene expression signature, Stem.Sig, that successfully predicts response to immunotherapy across multiple cancer types.
We found that Stem.Sig genes were enriched in the following biological functions: hypoxia, glycolysis, ubiquitination, nucleotide excision repair, EPH-ephrin signaling, and WNT signaling. WNT signaling is the key pathway that drives self-renewal of CSCs and maintains cancer stemness . Hypoxia causes an increase in transcription factors (e.g., OCT4, SOX2, c-myc, and Nanog) which contribute to the sustenance of CSCs . Anaerobic glycolysis is the distinct metabolic hallmark of stem cells . Ubiquitination-mediated transcriptional regulatory network is essential in the maintenance of the stemness and pluripotency of stem cells . Nucleotide excision repair (NER) is a major DNA repair pathway, which preserves genome integrity of cancer stem cells as to overcome stressful conditions . Activity of EPH-ephrin signaling, as the largest family of receptor tyrosine kinases, is found enhanced in CSCs . In our previous study, nonsynonymous somatic mutations of EPHA3 and EPHA7 was found associated with improved ICI efficacy . It is reasonable that elevated EPH-ephrin signaling may contribute to the immunosuppressive features of CSCs. Furthermore, we evaluated the correlation between Stem.Sig and twelve previously identified stemness signatures . As expected, Stem.Sig was found positively associated with these stemness signatures across different cancer types (Additional file 2: Fig. S3). Our results were in line with previous studies and suggested that Stem.Sig encompasses genes that robustly and specifically correlate with cancer stemness.
TCGA pan-caner transcriptomic analysis revealed a consistently down-regulated expression of immune-related genes and reduced infiltration of immune cells in tumors with high Stem.Sig level across different cancer types. Interestingly, a negative association between B cells and Stem.Sig was also observed. B cells could favorably affect ICI response via tertiary lymphoid structure (TLS), and hence we analyzed the relationship between TLS and Stem.Sig . TLS scores were found inversely associated with Stem.Sig (Additional file 2: Fig. S4). Further analysis also revealed an up-regulation of some immune-relevant biological functions, including metabolism, DNA repair, and MYC signaling. Acquisition of hypermetabolic phenotype is an evolving mechanism that mediates immune evasion . Enhanced DNA-repair capacity prepared malignant cells for unfriendly environments . Increased MYC signaling suppresses immune response by elevating expression of PD-L1 and CD47 . Tumors with high Stem.Sig presented with substantially immunosuppressive features, which corroborate the predictive value of Stem.Sig.
Also, we observed a positive correlation between Stem.Sig and both TMB and ITH, which is similar to the results of Miranda et al. . It is noteworthy that high TMB is associated with high stemness. Although TMB is a well-recognized ICI biomarker, there is still a significant number of patients with high TMB fail to response to ICI . Our stratified analysis revealed a significantly negative correlation between Stem.Sig and anti-tumor immunity in both low TMB and high TMB tumors. Cancer stemness can be a reasonable explanation of the immune resistance of high TMB tumors, which further stressed the importance of Stem.Sig as a predictive ICI biomarker.
Stem.Sig is a novel biomarker that is capable of predicting ICI response effectively and distinguishing patients with survival benefits successfully. We further compared Stem.Sig with other state-of-the-art signatures, including six pan-cancer signatures [76,77,78,79,80] and seven melanoma-specific signatures [7, 14, 81,82,83,84,85]. Stem.Sig outperformed pan-cancer signatures with better generalization and achieved an overall favorable performance in different cohorts across multiple cancer types. Compared with melanoma-specific signatures, Stem.Sig ranked top 3 and achieved a competitive AUC of 0.76.
Biomarker research is not only for improving patient selection but also for combination strategies that can overcome immune resistance. Considering such a robust link between Stem.Sig and ICI outcomes, we used CRISPR datasets to explore potential drug targets from Stem.Sig. We ranked genes based on their relevance to immune response and harvested the most immune-resistant Stem.Sig genes. For example, BECN1 is among the top-ranked Stem.Sig genes to render the TIME resistant to ICI. BECN1 plays a central role in autophagy, which is essential for self-renewal of CSCs and the maintenance of cancer stemness . Targeting BECN1 can induce expression of CCL5, promote infiltration of NK cells, and thus improve antitumor immune response . Top-ranked Stem.Sig genes, such as EMC3, BECN1, VPS35, and PCBP2, showed improved immune response after knockout in melanoma, renal carcinoma, breast carcinoma, and colon adenocarcinoma from multiple CRISPR datasets. These stemness-associated genes could be potential therapeutic targets for various cancer types. Further research of these top-ranked Stem.Sig genes would help to develop a combined strategy of immunotherapy.
Our study has some limitations. First, there were only treatment naïve patients and non-responders from GSE115978 . Comparison of the cancer stemness was conducted between non-responders and treatment naïve patients. Considering the average response rate of melanoma is 30–40%, a considerable proportion of treatment naïve patients would probably not response to ICI. Theoretically, the difference between TN and NR is smaller than that of R and NR, since TN is a mixture of NR and R. However, a significant difference of stemness level still existed between NR and TN in this study, which indicates an even greater gap between NR and R. And this was confirmed by analysis of another scRNA-Seq ICI cohort, GSE123813 . Secondly, some clinical annotation data (e.g., sex/age/tumor stage/TMB/ITH) was unavailable in some RNA-Seq ICI studies for multivariate cox regression analysis of overall survival. Thirdly, the 10 RNA-Seq ICI cohorts adopted in our studies only cover five cancer types (GC, SKCM, RCC, UC, and GBM). The consistent negative association between Stem.Sig and anti-tumor immunity across 30 cancer types can compensate this to some degree. Still, the predictive value of Stem.Sig in a pan-caner setting needs to be verified by future prospective ICI trials.
We provided the first solid clinical evidence that cancer stemness was associated with immunotherapy resistance. Using pan-cancer analysis of single-cell transcriptomic data, we developed a gene expression signature, Stem.Sig, which outweighs other well-established signatures in predicting ICI outcomes across multiple cohorts. Further exploration of Stem.Sig also revealed some potential therapeutic targets. Our study demonstrates a promising solution for patient selection in immunotherapy and sheds light on tackling ICI resistance through targeting cancer stemness to boost anti-tumor immunity.
Availability of data and materials
All data used in this study are publicly available as described in the Method section. The web links or unique identifiers for public cohorts/datasets are described in the paper. Source codes and all supplementary data used to generate the results were deposited on figshare (https://doi.org/10.6084/m9.figshare.17654633) .
Area under the ROC curve
Basal cell carcinoma
Cancer stem cells
Large B cell lymphoma
False discovery rate
Gene set variation analysis
Head and neck cancer
High Stem.Sig and high TMB
High Stem.Sig and low TMB
Immune checkpoint inhibitor
Acute myeloid leukemia
Liver hepatocellular carcinoma
Low Stem.Sig and high TMB
Low Stem.Sig and low TMB
Merkel cell carcinoma
Nucleotide excision repair
Non-small cell lung cancer
Over Representation Analysis
Objective response rate
Ovarian serous cystadenocarcinoma
Renal cell carcinoma
Response Evaluation Criteria in Solid Tumors
The receiver operating characteristic
Single-cell RNA sequencing
Skin cutaneous melanoma
The Cancer Genome Atlas
Total mutational burden
Whole exome sequencing
Wang Y, Wang M, Wu H, Xu R. Advancing to the era of cancer immunotherapy. Cancer Commun. 2021:cac2.12178. https://doi.org/10.1002/cac2.12178.
Sharma P, Siddiqui BA, Anandhan S, Yadav SS, Subudhi SK, Gao J, et al. The Next Decade of Immune Checkpoint Therapy. Cancer Discov. 2021;11:838–57. https://doi.org/10.1158/2159-8290.CD-20-1680.
Davoli T, Uno H, Wooten EC, Elledge SJ. Tumor aneuploidy correlates with markers of immune evasion and with reduced response to immunotherapy. Science. 2017;355:eaaf8399. https://doi.org/10.1126/science.aaf8399.
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:510–25. https://doi.org/10.1158/2159-8290.CD-18-0957.
Ott PA, Bang Y-J, Piha-Paul SA, Razak ARA, Bennouna J, Soria J-C, et al. T-Cell–Inflamed Gene-Expression Profile, Programmed Death Ligand 1 Expression, and Tumor Mutational Burden Predict Efficacy in Patients Treated With Pembrolizumab Across 20 Cancers: KEYNOTE-028. J Clin Oncol. 2019;37:318–27. https://doi.org/10.1200/JCO.2018.78.2276.
Zhang Z, Wu H, Lin W, Wang Z, Yang L, Zeng Z, et al. EPHA7 mutation as a predictive biomarker for immune checkpoint inhibitors in multiple cancers. BMC Med. 2021;19:26. https://doi.org/10.1186/s12916-020-01899-x.
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:1545–9. https://doi.org/10.1038/s41591-018-0157-9.
Fu J, Li K, Zhang W, Wan C, Zhang J, Jiang P, et al. Large-scale public data reuse to model immunotherapy response and resistance. Genome Med. 2020;12:21. https://doi.org/10.1186/s13073-020-0721-z.
Hwang B, Lee JH, Bang D. Single-cell RNA sequencing technologies and bioinformatics pipelines. Exp Mol Med. 2018;50:1–14. https://doi.org/10.1038/s12276-018-0071-8.
Chen P, Hsu W, Han J, Xia Y, Depinho RA. Cancer stemness meets immunity: from mechanism to therapy. CellReports. 2021;34:108597. https://doi.org/10.1016/j.celrep.2020.108597.
Bayik D, Lathia JD. Cancer stem cell-immune cell crosstalk in tumour progression. Nat Rev Cancer. 2021. https://doi.org/10.1038/s41568-021-00366-w.
Miranda A, Hamilton PT, Zhang AW, Pattnaik S, Becht E, Mezheyeuski A, et al. Cancer stemness, intratumoral heterogeneity, and immune response across cancers. Proc Natl Acad Sci. 2019;116:9020–9. https://doi.org/10.1073/pnas.1818210116.
Gulati GS, Sikandar SS, Wesche DJ, Manjunath A, Bharadwaj A, Berger MJ, et al. Single-cell transcriptional diversity is a hallmark of developmental potential. Science. 2020;367:405–11. https://doi.org/10.1126/science.aax0249.
Jerby-Arnon L, Shah P, Cuoco MS, Rodman C, Su M-J, Melms JC, et al. A Cancer Cell Program Promotes T Cell Exclusion and Resistance to Checkpoint Blockade. Cell. 2018;175:984–97.e24. https://doi.org/10.1016/j.cell.2018.09.006.
Yost KE, Satpathy AT, Wells DK, Qi Y, Wang C, Kageyama R, et al. Clonal replacement of tumor-specific T cells following PD-1 blockade. Nat Med. 2019;25:1251–9. https://doi.org/10.1038/s41591-019-0522-3.
Wang L, Dai J, Han R-R, Dong L, Feng D, Zhu G, et al. Single-cell map of diverse immune phenotypes in the metastatic brain tumor microenvironment of non small cell lung cancer. bioRxiv. 2019. https://doi.org/10.1101/2019.12.30.890517.
Venteicher AS, Tirosh I, Hebert C, Yizhak K, Neftel C, Filbin MG, et al. Decoupling genetics, lineages, and microenvironment in IDH-mutant gliomas by single-cell RNA-seq. Science. 2017;355. https://doi.org/10.1126/science.aai8478.
Tirosh I, Venteicher AS, Hebert C, Escalante LE, Patel AP, Yizhak K, et al. Single-cell RNA-seq supports a developmental hierarchy in human oligodendroglioma. Nature. 2016;539:309–13. https://doi.org/10.1038/nature20123.
Tirosh I, Izar B, Prakadan SM, Wadsworth MH, Treacy D, Trombetta JJ, et al. Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq. Science. 2016;352:189–96. https://doi.org/10.1126/science.aad0501.
Song Q, Hawkins GA, Wudel L, Chou PC, Forbes E, Pullikuth AK, et al. Dissecting intratumoral myeloid cell plasticity by single cell RNA-seq. Cancer Med. 2019;8:3072–85. https://doi.org/10.1002/cam4.2113.
Shih AJ, Menzin A, Whyte J, Lovecchio J, Liew A, Khalili H, et al. Correction: Identification of grade and origin specific cell populations in serous epithelial ovarian cancer by single cell RNA-seq. PLoS One. 2018;13:1–17. https://doi.org/10.1371/journal.pone.0208778.
Rao M, Oh K, Moffitt R, Thompson P, Li J, Liu J, et al. Comparative single-cell RNA sequencing (scRNA-seq) reveals liver metastasis–specific targets in a patient with small intestinal neuroendocrine cancer. Mol Case Stud. 2020;6:a004978. https://doi.org/10.1101/mcs.a004978.
Puram SV, Tirosh I, Parikh AS, Patel AP, Yizhak K, Gillespie S, et al. Single-Cell Transcriptomic Analysis of Primary and Metastatic Tumor Ecosystems in Head and Neck Cancer. Cell. 2017;171:1611–24.e24. https://doi.org/10.1016/j.cell.2017.10.044.
Peng J, Sun BF, Chen CY, Zhou JY, Chen YS, Chen H, et al. Single-cell RNA-seq highlights intra-tumoral heterogeneity and malignant progression in pancreatic ductal adenocarcinoma. Cell Res. 2019;29:725–38. https://doi.org/10.1038/s41422-019-0195-y.
Paulson KG, Voillet V, McAfee MS, Hunter DS, Wagener FD, Perdicchio M, et al. Acquired cancer resistance to combination immunotherapy from transcriptional loss of class I HLA. Nat Commun. 2018;9. https://doi.org/10.1038/s41467-018-06300-3.
Zilionis R, Engblom C, Pfirschke C, Savova V, Zemmour D, Saatcioglu HD, et al. Single-Cell Transcriptomics of Human and Mouse Lung Cancers Reveals Conserved Myeloid Populations across Individuals and Species. Immunity. 2019;50:1317–34.e10. https://doi.org/10.1016/j.immuni.2019.03.009.
Neftel C, Laffy J, Filbin MG, Hara T, Shore ME, Rahme GJ, et al. An Integrative Model of Cellular States, Plasticity, and Genetics for Glioblastoma. Cell. 2019;178:835–49.e21. https://doi.org/10.1016/j.cell.2019.06.024.
Moncada R, Barkley D, Wagner F, Chiodin M, Devlin JC, Baron M, et al. Integrating microarray-based spatial transcriptomics and single-cell RNA-seq reveals tissue architecture in pancreatic ductal adenocarcinomas. Nat Biotechnol. 2020;38:333–42. https://doi.org/10.1038/s41587-019-0392-8.
Ma L, Hernandez MO, Zhao Y, Mehta M, Tran B, Kelly M, et al. Tumor cell biodiversity drives microenvironmental reprogramming in liver cancer. Cancer Cell. 2019;36:418–430.e6. https://doi.org/10.1016/j.ccell.2019.08.007.
Ledergor G, Weiner A, Zada M, Wang S-Y, Cohen YC, Gatt ME, et al. Single cell dissection of plasma cell heterogeneity in symptomatic and asymptomatic myeloma. Nat Med. 2018;24:1867–76. https://doi.org/10.1038/s41591-018-0269-2.
Lambrechts D, Wauters E, Boeckx B, Aibar S, Nittner D, Burton O, et al. Phenotype molding of stromal cells in the lung tumor microenvironment. Nat Med. 2018;24:1277–89. https://doi.org/10.1038/s41591-018-0096-5.
Kim C, Gao R, Sei E, Brandt R, Hartman J, Hatschek T, et al. Chemoresistance evolution in triple-negative breast cancer delineated by single-cell sequencing. Cell. 2018;173:879–93.e13. https://doi.org/10.1016/j.cell.2018.03.041.
Hovestadt V, Smith KS, Bihannic L, Filbin MG, Shaw MKL, Baumgartner A, et al. Resolving medulloblastoma cellular architecture by single-cell genomics. Nature. 2019;572:74–9. https://doi.org/10.1038/s41586-019-1434-6.
Filbin MG, Tirosh I, Hovestadt V, Shaw ML, Escalante LE, Mathewson ND, et al. Developmental and oncogenic programs in H3K27M gliomas dissected by single-cell RNA-seq. Science. 2018;360:331–5. https://doi.org/10.1126/science.aao4750.
Durante MA, Rodriguez DA, Kurtenbach S, Kuznetsov JN, Sanchez MI, Decatur CL, et al. Single-cell analysis reveals new evolutionary complexity in uveal melanoma. Nat Commun. 2020;11. https://doi.org/10.1038/s41467-019-14256-1.
Darmanis S, Sloan SA, Croote D, Mignardi M, Chernikova S, Samghababi P, et al. Single-Cell RNA-Seq analysis of infiltrating neoplastic cells at the migrating front of human glioblastoma. Cell Rep. 2017;21:1399–410. https://doi.org/10.1016/j.celrep.2017.10.030.
Zhao W, Dovas A, Spinazzi EF, Levitin HM, Banu MA, Upadhyayula P, et al. Deconvolution of cell type-specific drug responses in human tumor tissue with single-cell RNA-seq. Genome Med. 2021;13. https://doi.org/10.1186/s13073-021-00894-y.
Zhang P, Yang M, Zhang Y, Xiao S, Lai X, Tan A, et al. Dissecting the single-cell transcriptome network underlying gastric premalignant lesions and early gastric cancer. Cell Rep. 2019;27:1934–47.e5. https://doi.org/10.1016/j.celrep.2019.04.052.
Zhang L, Li Z, Skrzypczynska KM, Fang Q, Zhang W, O’Brien SA, et al. Single-cell analyses inform mechanisms of myeloid-targeted therapies in colon cancer. Cell. 2020;181:442–59.e29. https://doi.org/10.1016/j.cell.2020.03.048.
Yuan J, Levitin HM, Frattini V, Bush EC, Boyett DM, Samanamud J, et al. Single-cell transcriptome analysis of lineage diversity in high-grade glioma. Genome Med. 2018;10:1–15. https://doi.org/10.1186/s13073-018-0567-9.
Wang R, Sharma R, Shen X, Laughney AM, Funato K, Clark PJ, et al. Adult human glioblastomas harbor radial glia-like cells. Stem Cell Rep. 2020;14:338–50. https://doi.org/10.1016/j.stemcr.2020.01.007.
Wang L, Catalan F, Shamardani K, Babikir H, Diaz A. Ensemble learning for classifying single-cell data and projection across reference atlases. Bioinformatics. 2020;36:3585–7. https://doi.org/10.1093/bioinformatics/btaa137.
Wang L, Babikir H, Müller S, Yagnik G, Shamardani K, Catalan F, et al. The phenotypes of proliferating glioblastoma cells reside on a single axis of variation. Cancer Discov. 2019;9:1708–19. https://doi.org/10.1158/2159-8290.CD-19-0329.
Goldman MJ, Craft B, Hastie M, Repečka K, McDade F, Kamath A, et al. Visualizing and interpreting cancer genomics data via the Xena platform. Nat Biotechnol. 2020;38:675–8. https://doi.org/10.1038/s41587-020-0546-8.
Lawson KA, Sousa CM, Zhang X, Kim E, Akthar R, Caumanns JJ, et al. Functional genomic landscape of cancer-intrinsic evasion of killing by T cells. Nature. 2020;586:120–6. https://doi.org/10.1038/s41586-020-2746-2.
Kearney CJ, Vervoort SJ, Hogg SJ, Ramsbottom KM, Freeman AJ, Lalaoui N, et al. Tumor immune evasion arises through loss of TNF sensitivity. Sci Immunol. 2018;3:eaar3451. https://doi.org/10.1126/sciimmunol.aar3451.
Freeman AJ, Vervoort SJ, Ramsbottom KM, Kelly MJ, Michie J, Pijpers L, et al. Natural killer cells suppress T cell-associated tumor immune evasion. Cell Rep. 2019;28:2784–94.e5. https://doi.org/10.1016/j.celrep.2019.08.017.
Vredevoogd DW, Kuilman T, Ligtenberg MA, Boshuizen J, Stecker KE, de Bruijn B, et al. augmenting immunotherapy impact by lowering tumor TNF cytotoxicity threshold. Cell. 2019;178:585–99.e15. https://doi.org/10.1016/j.cell.2019.06.014.
Patel SJ, Sanjana NE, Kishton RJ, Eidizadeh A, Vodnala SK, Cam M, et al. Identification of essential genes for cancer immunotherapy. Nature. 2017;548:537–42. https://doi.org/10.1038/nature23477.
Pan D, Kobayashi A, Jiang P, Ferrari de Andrade L, Tay RE, Luoma AM, et al. A major chromatin regulator determines resistance of tumor cells to T cell–mediated killing. Science. 2018;359:770–5. https://doi.org/10.1126/science.aao1710.
Manguso RT, Pope HW, Zimmer MD, Brown FD, Yates KB, Miller BC, et al. In vivo CRISPR screening identifies Ptpn2 as a cancer immunotherapy target. Nature. 2017;547:413–8. https://doi.org/10.1038/nature23270.
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:909–18. https://doi.org/10.1038/s41591-020-0839-y.
Mariathasan S, Turley SJ, Nickles D, Castiglioni A, Yuen K, Wang Y, et al. TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature. 2018;554:544–8. https://doi.org/10.1038/nature25501.
Liu D, Schilling B, Liu D, Sucker A, Livingstone E, Jerby-amon L, et al. Integrative molecular and clinical modeling of clinical outcomes to PD1 blockade in patients with metastatic melanoma. Nat Med. 2019;25. https://doi.org/10.1038/s41591-019-0654-5.
Gide TN, Quek C, Menzies AM, Tasker AT, Shang P, Holst J, et al. Distinct immune cell populations define response to Anti-PD-1 monotherapy and anti-PD-1/Anti-CTLA-4 combined therapy. Cancer Cell. 2019;35:238–55.e6. https://doi.org/10.1016/j.ccell.2019.01.003.
Riaz N, Havel JJ, Makarov V, Desrichard A, Urba WJ, Sims JS, et al. Tumor and microenvironment evolution during immunotherapy with nivolumab. Cell. 2017;171:934–49.e16. https://doi.org/10.1016/j.cell.2017.09.028.
Zhao J, Chen AX, Gartrell RD, Silverman AM, Aparicio L, Chu T, et al. Immune and genomic correlates of response to anti-PD-1 immunotherapy in glioblastoma. Nat Med. 2019;25:462–9. https://doi.org/10.1038/s41591-019-0349-y.
Snyder A, Nathanson T, Funt SA, Ahuja A, Buros Novik J, Hellmann MD, et al. Contribution of systemic and somatic factors to clinical response and resistance to PD-L1 blockade in urothelial cancer: an exploratory multi-omic analysis. PLoS Med. 2017;14:e1002309. https://doi.org/10.1371/journal.pmed.1002309 Minna JD, editor.
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:35–44. https://doi.org/10.1016/j.cell.2016.02.065.
Van Allen EM, Miao D, Schilling B, Shukla SA, Blank C, Zimmer L, et al. Erratum for the Report “Genomic correlates of response to CTLA-4 blockade in metastatic melanoma”. Science. 2016;352:aaf8264. https://doi.org/10.1126/science.aaf8264 by E. M. Van Allen, D. Miao, B. Schilling, S. A. Shukla, C. Blank, L. Zimmer, A. Sucker, U. Hillen, M. H. Geukes Foppen, S. M. Goldinger, J. Utikal, J. C. Ha.
Kim ST, Cristescu R, Bass AJ, Kim K-M, Odegaard JI, Kim K, et al. Comprehensive molecular characterization of clinical responses to PD-1 inhibition in metastatic gastric cancer. Nat Med. 2018;24:1449–58. https://doi.org/10.1038/s41591-018-0101-z.
Sun D, Wang J, Han Y, Dong X, Ge J, Zheng R, et al. TISCH: a comprehensive web resource enabling interactive single-cell transcriptome visualization of tumor microenvironment. Nucleic Acids Res. 2021;49:D1420–30. https://doi.org/10.1093/nar/gkaa1020.
Chen Y-T, Shen J-Y, Chen D-P, Wu C-F, Guo R, Zhang P-P, et al. Identification of cross-talk between m6A and 5mC regulators associated with onco-immunogenic features and prognosis across 33 cancer types. J Hematol Oncol. 2020;13:22. https://doi.org/10.1186/s13045-020-00854-w.
Gao J, Aksoy BA, Dogrusoz U, Dresdner G, Gross B, Sumer SO, et al. Integrative Analysis of Complex Cancer Genomics and Clinical Profiles Using the cBioPortal. Sci Signal. 2013;6:pl1. https://doi.org/10.1126/scisignal.2004088.
Cerami E, Gao J, Dogrusoz U, Gross BE, Sumer SO, Aksoy BA, et al. The cBio Cancer Genomics Portal: An Open Platform for Exploring Multidimensional Cancer Genomics Data: Figure 1. Cancer Discov. 2012;2:401–4. https://doi.org/10.1158/2159-8290.CD-12-0095.
Thorsson V, Gibbs DL, Brown SD, Wolf D, Bortone DS, Ou Yang T-H, et al. The Immune Landscape of Cancer. Immunity. 2018;48:812–30.e14. https://doi.org/10.1016/j.immuni.2018.03.023.
Boyle EI, Weng S, Gollub J, Jin H, Botstein D, Cherry JM, et al. GO::TermFinder--open source software for accessing Gene Ontology information and finding significantly enriched Gene Ontology terms associated with a list of genes. Bioinformatics. 2004;20:3710–5. https://doi.org/10.1093/bioinformatics/bth456.
Yu G, Wang L-G, Han Y, He Q-Y. clusterProfiler: an R Package for Comparing Biological Themes Among Gene Clusters. Omi A J Integr Biol. 2012;16:284–7. https://doi.org/10.1089/omi.2011.0118.
Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-Seq data. BMC Bioinformatics. 2013;14:7. https://doi.org/10.1186/1471-2105-14-7.
Becht E, Giraldo NA, Lacroix L, Buttard B, Elarouci N, Petitprez F, et al. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol. 2016;17:1–20. https://doi.org/10.1186/s13059-016-1070-5.
Eisenhauer EA, Therasse P, Bogaerts J, Schwartz LH, Sargent D, Ford R, et al. New response evaluation criteria in solid tumours: Revised RECIST guideline (version 1.1). Eur J Cancer. 2009;45:228–47. https://doi.org/10.1016/j.ejca.2008.10.026.
Johnson WE, Li C, Rabinovic A. Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics. 2007;8:118–27. https://doi.org/10.1093/biostatistics/kxj037.
Vougas K, Sakellaropoulos T, Kotsinas A, Foukas G-RP, Ntargaras A, Koinis F, et al. Machine learning and data mining frameworks for predicting drug response in cancer: An overview and a novel in silico screening process based on association rule mining. Pharmacol Ther. 2019;203:107395. https://doi.org/10.1016/j.pharmthera.2019.107395.
Budczies J, Kosztyla D, von Törne C, Stenzinger A, Drab-Esfahani S, Dietel M, et al. cancerclass: An R Package for development and validation of diagnostic tests from high-dimensional molecular data. J Stat Softw. 2014;59:1–8. https://doi.org/10.18637/jss.v059.i01.
Sammut S, Crispin-Ortuzar M, Chin S, Provenzano E, Bardwell HA, Ma W, et al. Multi-omic machine learning predictor of breast cancer therapy response. Nature. 2021. https://doi.org/10.1038/s41586-021-04278-5.
Ayers M, Lunceford J, Nebozhyn M, Murphy E, Loboda A, Kaufman DR, et al. IFN-γ–related mRNA profile predicts clinical response to PD-1 blockade. J Clin Invest. 2017;127:2930–40. https://doi.org/10.1172/JCI91190.
Topalian SL, Hodi FS, Brahmer JR, Gettinger SN, Smith DC, McDermott DF, et al. Safety, activity, and immune correlates of Anti–PD-1 antibody in cancer. N Engl J Med. 2012;366:2443–54. https://doi.org/10.1056/nejmoa1200690.
Dominguez CX, Müller 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:232–53. https://doi.org/10.1158/2159-8290.CD-19-0644.
Ju M, Bi J, Wei Q, Jiang L, Guan Q, Zhang M, et al. Pan-cancer analysis of NLRP3 inflammasome with potential implications in prognosis and immunotherapy in human cancer. Brief Bioinform. 2020;00:1–16. https://doi.org/10.1093/bib/bbaa345.
Rooney MS, Shukla SA, Wu CJ, Getz G, Hacohen N. Molecular and genetic properties of tumors associated with local immune cytolytic activity. Cell. 2015;160:48–61. https://doi.org/10.1016/j.cell.2014.12.033.
Shukla SA, Bachireddy P, Schilling B, Galonska C, Zhan Q, Bango C, et al. Cancer-germline antigen expression discriminates clinical outcome to CTLA-4 Blockade. Cell. 2018;173:624–33.e8. https://doi.org/10.1016/j.cell.2018.03.026.
Hugo W, Zaretsky JM, Sun L, Johnson DB, Ribas A, Lo RS, et al. Genomic and transcriptomic features of response to Anti-PD-1 therapy in metastatic melanoma article genomic and transcriptomic features of response to Anti-PD-1 therapy in metastatic melanoma. Cell. 2016:1–10. https://doi.org/10.1016/j.cell.2016.02.065.
Xiong D, Wang Y, You M. A gene expression signature of TREM2hi macrophages and γδ T cells predicts immunotherapy response. Nat Commun. 2020;11:1–12. https://doi.org/10.1038/s41467-020-18546-x.
Cui C, Xu C, Yang W, Chi Z, Sheng X, Si L, et al. Ratio of the interferon-γ signature to the immunosuppression signature predicts anti-PD-1 therapy response in melanoma. npj. Genomic Med. 2021;6:7. https://doi.org/10.1038/s41525-021-00169-w.
Yan M, Hu J, Ping Y, Xu L, Liao G, Jiang Z, et al. Single-cell transcriptomic analysis reveals a tumor-reactive T Cell signature associated with clinical outcome and immunotherapy response in melanoma. Front Immunol. 2021;12:1–14. https://doi.org/10.3389/fimmu.2021.758288.
Fluss R, Faraggi D, Reiser B. Estimation of the youden index and its associated cutoff point. Biom J. 2005;47:458–72. https://doi.org/10.1002/bimj.200410135.
Akoglu H. User’s guide to correlation coefficients. Turkish J Emerg Med. 2018;18:91–3. https://doi.org/10.1016/j.tjem.2018.08.001.
Xiao Q, Wu J, Wang WJ, Chen S, Zheng Y, Yu X, et al. DKK2 imparts tumor immunity evasion through β-catenin-independent suppression of cytotoxic immune-cell activation. Nat Med. 2018;24:262–70. https://doi.org/10.1038/nm.4496.
Rosenbaum M, Gewies A, Pechloff K, Heuser C, Engleitner T, Gehring T, et al. Bcl10-controlled Malt1 paracaspase activity is key for the immune suppressive function of regulatory T cells. Nat Commun. 2019;10. https://doi.org/10.1038/s41467-019-10203-2.
Richter C, Mayhew D, Rennhack JP, So J, Stover EH, Hwang JH, et al. Genomic amplification and functional dependency of the gamma actin gene ACTG1 in uterine cancer. Int J Mol Sci. 2020;21:8690. https://doi.org/10.3390/ijms21228690.
Mandili G, Curcio C, Bulfamante S, Follia L, Ferrero G, Mazza E, et al. In pancreatic cancer, chemotherapy increases antitumor responses to tumor-associated antigens and potentiates DNA vaccination. J Immunother Cancer. 2020;8:e001071. https://doi.org/10.1136/jitc-2020-001071.
Qin G, Wang X, Ye S, Li Y, Chen M, Wang S, et al. NPM1 upregulates the transcription of PD-L1 and suppresses T cell activity in triple-negative breast cancer. Nat Commun. 2020;11. https://doi.org/10.1038/s41467-020-15364-z.
Maugeri-Saccà M, Bartucci M, De Maria R. DNA damage repair pathways in cancer stem cells. Mol Cancer Ther. 2012;11:1627–36. https://doi.org/10.1158/1535-7163.MCT-11-1040.
Jaiswal AR, Liu AJ, Pudakalakatti S, Dutta P, Jayaprakash P, Bartkowiak T, et al. Melanoma evolves complete immunotherapy resistance through the acquisition of a hypermetabolic phenotype. Cancer Immunol Res. 2020;8:1365–80. https://doi.org/10.1158/2326-6066.CIR-19-0005.
Casey SC, Tong L, Li Y, Do R, Walz S, Fitzgerald KN, et al. MYC regulates the antitumor immune response through CD47 and PD-L1. Science. 2016;352:227–31. https://doi.org/10.1126/science.aac9935.
Wang Z, Wang Y, Yang T, Xing H, Wang Y, Gao L, et al. Machine learning revealed stemness features and a novel stemness-based classification with appealing implications in discriminating the prognosis, immunotherapy and temozolomide responses of 906 glioblastoma patients. Brief Bioinform. 2021;00:1–20. https://doi.org/10.1093/bib/bbab032.
Maccalli C, Rasul KI, Elawad M, Ferrone S. The role of cancer stem cells in the modulation of anti-tumor immune responses. Semin Cancer Biol. 2018;53:189–200. https://doi.org/10.1016/j.semcancer.2018.09.006.
Clara JA, Monge C, Yang Y, Takebe N. Targeting signalling pathways and the immune microenvironment of cancer stem cells — a clinical update. Nat Rev Clin Oncol. 2020;17:204–32. https://doi.org/10.1038/s41571-019-0293-2.
Abou Khouzam R, Goutham HV, Zaarour RF, Chamseddine AN, Francis A, Buart S, et al. Integrating tumor hypoxic stress in novel and more adaptable strategies for cancer immunotherapy. Semin Cancer Biol. 2020;65:140–54. https://doi.org/10.1016/j.semcancer.2020.01.003.
Wei P, Dove KK, Bensard C, Schell JC, Rutter J. The force is strong with this one: metabolism (Over)powers stem cell fate. Trends Cell Biol. 2018;28:551–9. https://doi.org/10.1016/j.tcb.2018.02.007.
Deng L, Meng T, Chen L, Wei W, Wang P. The role of ubiquitination in tumorigenesis and targeted drug discovery. Signal Transduct Target Ther. 2020;5. https://doi.org/10.1038/s41392-020-0107-0.
Chen J, Song W, Amato K. Eph receptor tyrosine kinases in cancer stem cells. Cytokine Growth Factor Rev. 2015;26:1–6. https://doi.org/10.1016/j.cytogfr.2014.05.001.
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:549–55. https://doi.org/10.1038/s41586-019-1922-8.
Ready N, Hellmann MD, Awad MM, Otterson GA, Gutierrez M, Gainor JF, et al. First-line nivolumab plus ipilimumab in advanced non–small-cell lung cancer (CheckMate 568): outcomes by programmed death ligand 1 and tumor mutational burden as biomarkers. J Clin Oncol. 2019;37:992–1000. https://doi.org/10.1200/JCO.18.01042.
Pan H, Cai N, Li M, Liu GH, Izpisua Belmonte JC. Autophagic control of cell “stemness”. EMBO Mol Med. 2013;5:327–31. https://doi.org/10.1002/emmm.201201999.
Mgrditchian T, Arakelian T, Paggetti J, Noman MZ, Viry E, Moussay E, et al. Targeting autophagy inhibits melanoma growth by enhancing NK cells infiltration in a CCL5-dependent manner. Proc Natl Acad Sci. 2017;114:E9271–9. https://doi.org/10.1073/pnas.1703921114.
Zhen Z, Zi-Xian W, Yan-Xing C, Hao-Xiang W, Ling Y, Qi Z, et al. Integrated analysis of single-cell and bulk RNA sequencing data reveals a pan-cancer stemness signature predicting immunotherapy response. Figshare. 2022. https://doi.org/10.6084/m9.figshare.17654633.
We would like to thank the staff members of the TCGA Research Network, the cBioPortal, the UCSC Xena data portal, and the TISCH data portal; as well as all the authors for making their valuable research data public.
This work was supported by grants from the Science and Technology Program of Guangdong (2019B020227002); Science and Technology Program of Guangzhou (201803040019, 20104020228, 202002030208); CAMS Innovation Fund for Medical Sciences (CIFMS) (2019-I2M-5-036); National Natural Science Foundation of China (81930065, 81871985, 82003269, 82073377, and 81772587); Natural Science Foundation of Guangdong Province (2014A030312015, 2019A1515011109, 2021A1515012439); China Postdoctoral Science Foundation (2021M693651).
Ethics approval and consent to participate
All datasets used in this study have been previously published.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Characteristics of ICIs scRNA cohorts. Table S2. List of scRNA datasets applied to develop Stem.Sig. Table S3. List of immunotherapy cohorts used in this study. Table S4. List of CRISPR datasets. Table S5. List of predictive gene expression signatures for immunotherapy. Table S6. Gene list of Stem.Sig. Table S7. Comparison of AUC of previous signatures in testing cohort.
Boxplots depicting the correlation of immune cells infiltration with Stem.Sig and TMB. Figure S2. Subgroup analysis of testing set. Figure S3. Correlation plot of Stem.Sig and other stemness signatures. Figure S4. Association Stem.Sig and expression level of TLS-related genes.
About this article
Cite this article
Zhang, Z., Wang, ZX., Chen, YX. et al. Integrated analysis of single-cell and bulk RNA sequencing data reveals a pan-cancer stemness signature predicting immunotherapy response. Genome Med 14, 45 (2022). https://doi.org/10.1186/s13073-022-01050-w