Skip to main content
Download PDF

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) Cite this article

Abstract

Background

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.

Methods

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.

Results

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.

Conclusions

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.

Background

Immune checkpoint inhibitor (ICI) has ushered in a new era of cancer treatment and provided unprecedented clinical benefits for patients [1]. However, only a relatively small proportion of patients respond to it [2], 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 [9].

Cancer stem cells (CSCs) are self-renewal cells that promote tumor initiation, progression, and metastasis [10]. Mounting evidences revealed a prominent association between stemness and cancer immune evasion and resistance [11]. A previous study demonstrated that high stemness correlates with immune cell exclusion across 21 solid cancer types [12], 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 [13].

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) [44], 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.

Methods

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 [14] and GSE123813 [15], 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/) [62], 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) [44] 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 [63]. Total mutation burden (TMB) was retrieved from cBioPortal (https://www.cbioportal.org) [64, 65] and intratumor heterogeneity (ITH) data was from Thorsson et al. [66], 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 [59], Liu 2019 [54], Gide 2019 [55], Riaz 2017 [56], Van Allen 2015 [60]), 2 urothelial carcinoma (UC) cohorts (Mariathasan 2018 [53], Synder 2017 [58]), 1 glioblastoma multiform (GBM) cohort (Zhao 2019 [57]), 1 gastric cancer (GC) cohort (Kim 2018 [61]) and 1 renal cell carcinoma (RCC) cohort (Braun 2020 [52]). 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 [59] comprises 27 pre-treated tumor samples from 26 patients, while cohort Zhao 2019 [57] 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 [47], Kearney 2018 [46], Manguso 2017 [51], Pan 2018 [50], Patel 2017 [49], Vredevoogd 2019 [48], and Lawson 2020 [45]. The first six CRISPR studies have been previously curated by Fu et al. [8] In addition to Fu et al., we further collected another CRISPR cohort from Lawson et al. [45]. 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) [45]. 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 [13]. 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 [13]. 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 [13]. 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 [62].

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 [67]. The R package clusterProfiler v4.2.1 was used to perform ORA and visualize the results [68].

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 [69].

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. [66]. The R package MCP-counter v1.1.0 was utilized to estimate the abundance of tumor-infiltrating leukocytes [70].

Clinical outcomes

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 [71], except cohort Hugo 2016 [59], 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

Dataset

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 [52], Mariathasan 2018 UC [53], Liu 2019 SKCM [54], Gide 2019 SKCM [55], and Riaz 2017 SKCM [56]. We used ComBat method to remove the batch effect of different ICI RNA-Seq cohorts [72]. 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 [75]. 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 [76], T.cell.inflamed.Sig [76], PDL1.Sig [77], LRRC15.CAF.Sig [78], NLRP3.Sig [79], and Cytotoxic.Sig [80]) and seven melanoma-specific signatures (CRMA.Sig [81], IMPRES.Sig [7], IPRES.Sig [82], TcellExc.Sig [14], ImmmunCells.Sig [83], IMS.Sig [84], and TRS.Sig [85]). 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 [59] and Van Allen 2015 [60]). Codes and algorithms for the 13 aforementioned signatures were derived from their original studies, such as ssGSEA for NLRP3.Sig [79], cancerclass for ImmuneCell.Sig [83], overall expression for TcellExc.Sig [14], and so on. Details of these signatures and their corresponding algorithms can be found in Additional file 1: Table S5.

Statistical analysis

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 [74]. 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 [86]. 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.

Results

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 [14]. 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 [14]. 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 [15]. 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).

Fig. 1
figure 1

Identification and validation of a negative association between cancer cell stemness and ICI outcomes. A, C t-Distributed Stochastic Neighbor Embedding (tSNE) plot of malignant cells from SKCM or BCC. Top tSNE plots depicting the distribution of CytoTRACE scores among malignant cells. Dark-green indicates lower scores (low stemness) while dark-red indicates higher scores (high stemness). Bottom tSNE plots label the malignant cells by response phenotype. B, D raincloud plot of CytoTRACE scores by response phenotype (NR vs. TN) in SKCM cohort or by response phenotype (NR vs. R) in BCC cohort. The center of the box plot are median values, the bounds of the box are 25% and 75% quantiles (Wilcoxon test; *** P < 0.001). Abbreviation: NR, non-responders; R, responders; TN, treatment naïve patients.

Full size image

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 [14]. 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 [87].

Fig. 2
figure 2

Development and description of stemness signature. A Circos plot depicting the development of Stem.Sig. B Pathway enrichment analysis of genes in Stem.Sig. The bar plot showed the top 20 enriched Reactome pathways. The cnetplot presented the network of specific genes from these pathways. Colored points referred to the corresponding pathways. Abbreviation: CFTR, cystic fibrosis transmembrane conductance regulator; GG-NER, global genomic nucleotide excision repair; HIF, hypoxia-inducible factor; PCP, planar cell polarity; CE, convergent extension

Full size image

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 [66]. 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.

Fig. 3
figure 3

Analysis of the potential links between Stem.Sig and immune resistance using pan-cancer TCGA cohort. A Circos plot depicting the correlation between Stem.Sig and the expression level of immune-related genes across multiple cancer types. From inside to outside of the circos plot, the vertical axis with a black arrow indicated different cancer types, which were annotated by the y axis of plot B. B Heatmap depicting the correlation between Stem.Sig and the infiltration of immune cells across multiple cancer types. C Heatmap depicting the correlation between Stem.Sig and the Top 10 Hallmark pathways. D Correlation of median Stem.Sig and median TMB of each cancer type. E Correlation of median Stem.Sig and median ITH of each cancer type. GSVA scores were calculated to estimate the expression level of Stem.Sig for each sample

Full size image

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 [12]. 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).

Fig. 4
figure 4

Prediction of ICI outcomes using Stem.Sig. A Flow chart of training, validating, and testing the Stem.Sig model constructed using machine learning process. In the training set, we applied 10-time repeated 5-fold cross-validation for parameters tuning of different machine learning algorithms. In the validation set, Naïve Bayes algorithm with best AUC was kept as the final Stem.Sig model. (parameter: fL=0; adjust = 0.75; useKernel = TRUE). B Comparison of multiple ROC plot depicting the performance of different machine learning algorithms in the validation set. C ROC plot depicting the performance of the final Stem.Sig model in validation and testing cohort. D Kaplan-Meier curves comparing OS between High-risk and Low-risk patients in validation and testing set. “NR” and “R” predicted by the final Stem.Sig Model was defined as “High-risk” and “Low-risk” patients respectively. HR were calculated by Cox proportional hazards regression analysis. Abbreviation: TPR, true positive rate; FPR, false positive rate; AUC, area under the curve; HR, hazard ratio; CI, confidence intervals

Full size image

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 [76], T.cell.inflamed.Sig [76], PDL1.Sig [77], LRRC15.CAF.Sig [78], NLRP3.Sig [79], and Cytotoxic.Sig [80]), 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 [81], IMPRES.Sig [7], IPRES.Sig [82], TcellExc.Sig [14], ImmmunCells.Sig [83], IMS.Sig [84], and TRS.Sig [85]), 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.

Fig. 5
figure 5

Comparing AUC of Stem.Sig with other predictive gene signatures. A Circos plot depicting the performance of pan-cancer signatures in the testing set. The vertical axis indicated AUC values. Testing set comprises five different cohorts, including Hugo 2020 SKCM, Van Allen 2015 SKCM, Kim 2018 GC, Zhao 2019 GBM, Synder 2017 UC. B Heatmap comparing the predictive value of Stem.Sig and other pan-cancer signatures. Different signature rows were ordered by their AUC in the testing set. From top to bottom, Stem.Sig ranked first while Cytotoxic.Sig ranked last. C Bar plot depicting the AUC values of Stem.Sig and other melanoma-specific signatures in the SKCM cohort (Hugo 2016 + Van Allen 2015).

Full size image

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.

Fig. 6
figure 6

Exploration of potential treatment targets from Stem.Sig using CRISPR screening data. A Ranking of genes based on their knockout effects on anti-tumor immunity across 17 CRISPR datasets. Negative (positive) z scores indicated better (worse) immune response after knockout of a specific gene. Genes were ranked according to their mean z scores. Top-ranking genes were associated with immune resistance. Blank squares in the heatmap referred to missing values of gene data from the corresponding cohort. B Radar plot comparing the percentage of top-ranked genes for Stem.Sig and other predictive signatures. C Heatmap depicting z scores of 20 Stem.Sig genes in the 3% top-ranked genes across different CRISPR datasets

Full size image

Discussion

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 [10]. 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 [98]. Hypoxia causes an increase in transcription factors (e.g., OCT4, SOX2, c-myc, and Nanog) which contribute to the sustenance of CSCs [99]. Anaerobic glycolysis is the distinct metabolic hallmark of stem cells [100]. Ubiquitination-mediated transcriptional regulatory network is essential in the maintenance of the stemness and pluripotency of stem cells [101]. Nucleotide excision repair (NER) is a major DNA repair pathway, which preserves genome integrity of cancer stem cells as to overcome stressful conditions [93]. Activity of EPH-ephrin signaling, as the largest family of receptor tyrosine kinases, is found enhanced in CSCs [102]. In our previous study, nonsynonymous somatic mutations of EPHA3 and EPHA7 was found associated with improved ICI efficacy [6]. 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 [12]. 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 [103]. 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 [94]. Enhanced DNA-repair capacity prepared malignant cells for unfriendly environments [93]. Increased MYC signaling suppresses immune response by elevating expression of PD-L1 and CD47 [95]. 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. [12]. 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 [104]. 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 [105]. Targeting BECN1 can induce expression of CCL5, promote infiltration of NK cells, and thus improve antitumor immune response [106]. 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 [14]. 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 [15]. 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.

Conclusions

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) [107].

Abbreviations

AUC:

Area under the ROC curve

BCC:

Basal cell carcinoma

B-H:

Benjamini-Hochberg procedure

BRCA:

Breast cancer

CHOL:

Cholangiocarcinoma

CR:

Complete response

CRC:

Colorectal cancer

CSCs:

Cancer stem cells

DLBC:

Large B cell lymphoma

FDR:

False discovery rate

GSVA:

Gene set variation analysis

HNSC:

Head and neck cancer

HSHT:

High Stem.Sig and high TMB

HSLT:

High Stem.Sig and low TMB

ICI:

Immune checkpoint inhibitor

IRGs:

Immune-related genes

ITH:

Intra-tumoral heterogenicity

LAML:

Acute myeloid leukemia

LIHC:

Liver hepatocellular carcinoma

logFC:

Log-fold change

LSHT:

Low Stem.Sig and high TMB

LSLT:

Low Stem.Sig and low TMB

MB:

Medulloblastoma

MCC:

Merkel cell carcinoma

MM:

Multiple myeloma

NER:

Nucleotide excision repair

NET:

Neuroendocrine tumor

NSCLC:

Non-small cell lung cancer

NR:

Non-responders

ORA:

Over Representation Analysis

ORR:

Objective response rate

OS:

Overall survival

OV:

Ovarian serous cystadenocarcinoma

PAAD:

Pancreatic adenocarcinoma

PD:

Progressive disease

PR:

Partial response

R:

Responders

RCC:

Renal cell carcinoma

RECIST:

Response Evaluation Criteria in Solid Tumors

ROC:

The receiver operating characteristic

scRNA-Seq:

Single-cell RNA sequencing

SD:

Stable disease

SKCM:

Skin cutaneous melanoma

STAD:

Stomach adenocarcinoma

TCGA:

The Cancer Genome Atlas

THYM:

Thymoma

TILs:

Tumor-infiltrating leukocytes

TMB:

Total mutational burden

TN:

Treatment-naïve patients

UVM:

Uveal melanoma

WES:

Whole exome sequencing

References

  1. 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.

  2. 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.

    Article  CAS  PubMed  Google Scholar 

  3. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  4. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. 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.

    Article  PubMed  Google Scholar 

  6. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  9. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. 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.

    Article  CAS  Google Scholar 

  11. 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.

  12. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. 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.

  17. 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.

  18. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. 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.

    Article  Google Scholar 

  22. 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.

    Article  CAS  Google Scholar 

  23. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. 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.

  26. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. 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.

    Article  CAS  PubMed  Google Scholar 

  29. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. 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.

    Article  CAS  PubMed  Google Scholar 

  31. 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.

    Article  CAS  PubMed  Google Scholar 

  32. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. 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.

  36. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. 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.

  38. 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.

    Article  CAS  PubMed  Google Scholar 

  39. 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.

    Article  CAS  PubMed  Google Scholar 

  40. 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.

    Article  CAS  Google Scholar 

  41. 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.

    Article  CAS  Google Scholar 

  42. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. 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.

    Article  PubMed  Google Scholar 

  47. 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.

    Article  CAS  PubMed  Google Scholar 

  48. 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.

    Article  CAS  PubMed  Google Scholar 

  49. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. 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.

  55. 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.

    Article  CAS  PubMed  Google Scholar 

  56. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. 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.

    Article  Google Scholar 

  61. 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.

    Article  CAS  PubMed  Google Scholar 

  62. 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.

    Article  CAS  PubMed  Google Scholar 

  63. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. 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.

    Article  PubMed  Google Scholar 

  66. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  67. 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.

    Article  CAS  PubMed  Google Scholar 

  68. 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.

    Article  CAS  Google Scholar 

  69. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  70. 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.

    Article  CAS  Google Scholar 

  71. 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.

    Article  CAS  PubMed  Google Scholar 

  72. 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.

    Article  PubMed  Google Scholar 

  73. 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.

    Article  CAS  PubMed  Google Scholar 

  74. 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.

    Article  Google Scholar 

  75. 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.

  76. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  77. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  78. 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.

    Article  CAS  PubMed  Google Scholar 

  79. 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.

    Article  CAS  Google Scholar 

  80. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  81. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  82. 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.

  83. 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.

    Article  CAS  Google Scholar 

  84. 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.

    Article  CAS  Google Scholar 

  85. 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.

    Article  Google Scholar 

  86. 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.

    Article  PubMed  Google Scholar 

  87. 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.

    Article  Google Scholar 

  88. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  89. 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.

  90. 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.

    Article  CAS  PubMed Central  Google Scholar 

  91. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  92. 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.

  93. 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.

    Article  CAS  PubMed  Google Scholar 

  94. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  95. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  96. 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.

    Article  CAS  Google Scholar 

  97. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  98. 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.

    Article  PubMed  Google Scholar 

  99. 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.

    Article  CAS  PubMed  Google Scholar 

  100. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  101. 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.

  102. 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.

    Article  CAS  PubMed  Google Scholar 

  103. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  104. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  105. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  106. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  107. 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.

Download references

Acknowledgements

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.

Funding

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).

Author information

Author notes

  1. Zhen Zhang, Zi-Xian Wang, Yan-Xing Chen and Hao-Xiang Wu contributed equally to this work.

Authors and Affiliations

  1. Department of Medical Oncology, Sun Yat-sen University Cancer Center, State Key Laboratory of Oncology in South China, Collaborative Innovation Center for Cancer Medicine, Sun Yat-sen University, Guangzhou, 510060, P. R. China

    Zhen Zhang, Zi-Xian Wang, Yan-Xing Chen, Hao-Xiang Wu, Ling Yin, Qi Zhao, Hui-Yan Luo, Zhao-Lei Zeng, Miao-Zhen Qiu & Rui-Hua Xu

  2. Research Unit of Precision Diagnosis and Treatment for Gastrointestinal Cancer, Chinese Academy of Medical Sciences, Guangzhou, 510060, P. R. China

    Zhen Zhang, Zi-Xian Wang, Yan-Xing Chen, Hao-Xiang Wu, Ling Yin, Qi Zhao, Hui-Yan Luo, Zhao-Lei Zeng, Miao-Zhen Qiu & Rui-Hua Xu

  3. Laboratory of Artificial Intelligence and Data Science, Sun Yat-sen University Cancer Center, Guangzhou, 510060, P. R. China

    Zi-Xian Wang, Hui-Yan Luo & Rui-Hua Xu

Authors
  1. Zhen Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Zi-Xian Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Yan-Xing Chen
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Hao-Xiang Wu
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Ling Yin
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Qi Zhao
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Hui-Yan Luo
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Zhao-Lei Zeng
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Miao-Zhen Qiu
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Rui-Hua Xu
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

RHX and MZQ designed and supervised the study as joint senior authors; ZZ, ZXW, YXC, and HXW contributed to data collection and curation. ZZ and ZXW analyzed the datasets and interpreted the results. All authors interpreted the data, contributed to the writing and review of the manuscript, had full access to all the data in the study, and had final consent for publication. The authors read and approved the final manuscript.

Corresponding authors

Correspondence to Miao-Zhen Qiu or Rui-Hua Xu.

Ethics declarations

Ethics approval and consent to participate

All datasets used in this study have been previously published.

Consent for publication

Not applicable

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Table S1.

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.

Additional file 2: Figure S1.

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.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

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

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13073-022-01050-w

Keywords

  • Big data analysis
  • Single-cell sequencing
  • Immune checkpoint therapy
  • Stemness
  • Pan-cancer

Genome Medicine

ISSN: 1756-994X

Contact us