Summary of included CRISPR screens
In this study, a total of 17 CRISPR screens on the regulators of immune-mediated tumor elimination were included (the other five ICB-treated screens were described in the following section) (Table S1). To improve comprehension, a schematic diagram describing the process of CRISPR screen in a tumor/immune co-culture system is shown in Fig. 1A. Five cancers were selected as research objects, among which skin cancer (41.18%) was the most commonly studied cancer model, followed by breast cancer (23.53%) and colon cancer (23.53%) (Fig. 1B). The distribution of screens among different screening types, library types, organisms, and algorithms was also analyzed. We found that most of the included screens were carried out under in vitro conditions (76.47%), using a genome-scale library (88.24%) and mouse-derived cell line models (94.12%). Moreover, DrugZ (47.06%) and MAGeCK (35.29%) were the two most adopted algorithms for data analysis (Fig. S1). Notably, the included screens all used CRISPR knockout (CRISPRko) libraries, which were designed to investigate the relationship between the loss-of-function of genes and corresponding phenotypes.
Determination of antitumor immunity-related regulators
Based on data from these screens, we sought to find out potential regulator genes involved in antitumor immunity. Each screen could identify a plethora of enriched or depleted single guide RNAs (sgRNAs). Considering that no gain-of-function screens were included in this study, results from all the screens here could be interpreted in the same manner. Briefly, the knockout of genes that exert positive effects on anti-tumor immune response (such as genes associated with antigen presentation) resulted in resistance of corresponding tumor cells to T-cell mediated killing, thus enriching their sgRNAs (these genes were named as sensitizer genes hereafter). Conversely, the knockout of genes that promote tumor immune escape (such as CD274) can enhance the sensitivity of tumor cells to T cell-mediated killing, leading to the depletion of sgRNAs (these genes were named as resistor genes hereafter). The sensitizer and resistor genes in each screen were identified respectively. A tumor type-dependent result can be found in Table S2. All the results from different screens were then integrated and only those genes presented across two or more screens were selected for further analyses (Fig. 1C). This step yielded a preliminary list of 181 sensitizers and 427 resistors.
Given that the overall accuracy of screening results remained far from perfect, we intended to further narrow down this gene list to obtain more reliable candidate regulators. To this end, three well-proven tumor immunity-related gene signatures, namely ESTIMATE-based immune signature [19], cytolytic activity (CYT) signature [20], and MHC signature [21], were leveraged to discern sensitizer and resistor genes with potential functional significance. Specifically, we first calculated the corresponding scores of these three signatures based on the expression profiles from the TCGA Pan-Cancer cohort [22]. The functional status (a dichotomous variable, 0 = non-inactivation; 1 = inactivation) of each candidate sensitizer/resistor gene in each TCGA sample was determined using multi-omics TCGA data, which included expression, mutation, and copy number variation (CNV) data (see Materials and methods section). Given that all the candidate regulators were derived from the results of CRISPRko screens, we defined loss-of-function (inactivating) status as the main functional event. Then, the associations between the functional status of each sensitizer/resistor gene and the abundance of each signature could be determined by adopting a regression-based approach. After controlling for cancer type and adjusting for the multiple testing, associations with adj. P < 0.05 were deemed statistically significant. It should be noted that higher scores obtained for the three signatures indicated an enhanced anti-tumor immune response [23]. Therefore, according to the definition of sensitizer/resistor genes, sensitizers should be negatively associated with these signatures (lower scores upon inactivation), while resistors should be positively associated (higher scores upon inactivation). As indicated by the results, resistors exhibited an incline for a higher proportion of positive associations, albeit not significant, which could prove the rationality of this filtering approach to some extent (Fig. 1D). Functional sensitizer genes were defined as genes with significant negative associations (adj. P < 0.05 and coefficient < 0) with all three signatures (n = 65), while functional resistor genes were defined as genes with significant positive associations (adj. P < 0.05 and coefficient > 0) with all three signatures (n = 40) (Fig. 1E, Table S3).
Functional analysis of sensitizers and resistors
We conducted gene ontology (GO)-based functional similarity (FS) and annotation analysis to determine the functional characteristics of these sensitizer and resistor genes. The FS scores of each gene pair across all 105 genes were first computed and visualized on a heatmap (Fig. 2A). Besides, the distribution of FS scores of each sensitizer and resistor gene was also analyzed (Fig. S2A). Conceptually, a gene that had higher similarity with others was more likely to possess a central or significant role. In this regard, some critical transcription factor genes involved in tumor immunity, such as STAT1 and STAT2, were observed to have relatively high similarities with other genes. Intriguingly, the comparison between sensitizer and resistor genes suggested that resistors had significantly higher internal similarity (P = 0.0001, Fig. S2B). To characterize the biological processes related to these sensitizers and resistors, we next performed functional annotation based on the occurrence frequency of GO-BP terms from the Molecular Signatures Database (MSigDB) [24]. The results suggested that sensitizer genes tended to engage in immune-related processes, while resistor genes were mostly involved in some metabolic and biosynthetic processes (Fig. 2B, Table S4).
Dysregulation patterns of sensitizers and resistors
We also sought to explore the dysregulated expression patterns of these genes across different cancer types and immune subtypes. First, we focused on investigating the expression pattern of genes between tumor tissues and adjacent tissues. Differential expression analysis was conducted on 105 sensitizer and resistor genes across cancers in both tumor and normal tissues (n > 2). The results suggested that the differential patterns of the two gene types were inconsistent across different cancer types. Overall, in 57.9% of all cancer types, a numerically higher proportion of up-regulated resistor genes could be observed (Fig. S3). Next, the distribution of the inactivation events of sensitizer and resistor genes across different cancer types was analyzed (Fig. 2C). A more objective relationship was observed between resistors and sensitizers in each cancer type by calculating the normalized resistor-sensitizer ratio (Fig. S4A).
In addition, we also delineated the relationship between the dichotomous functional status of these genes and multi-class immune subtypes from a previous publication (C1-C6) [25]. To identify subtype-specific sensitizers or resistors, a logistic regression-based approach was adopted, and genes with adj. P < 0.05 and log2 odds ratio (OR) > 1.5 in certain immune subtypes were defined as the specific genes of this subtype. The analysis yielded a total of 37 subtype-specific genes (35 out of 37 were sensitizers), including two C1-specific genes, 10 C4-specific genes, and 25 C5-specific genes (Fig. 2D). The functions of these genes were then annotated using GO information from ImmPort (immport.org) [26]. Interestingly, we observed that the inactivation events of multiple C5-specific genes were related to the process of antigen processing and presentation (Fig. 2E). Then, the distribution of the inactivation events of sensitizers and resistors across six immune subtypes was analyzed (Fig. 2F), and the corresponding normalized sensitizer-resistor ratios were also calculated (Fig. S4B). Among the six subtypes, C5 had the highest average number of inactivation events of sensitizers (n = 9.24) as well as the highest normalized sensitizer-resistor ratio (r = 2.89), which was consistent with its property as an immunologically quiet subtype [25].
The functional characterization of sensitizers and resistors across cancers
Theoretically, the loss of sensitizer genes (such as B2M) should enable tumors to resist immune attack, while the loss of resistor genes (such as CD274) could augment the cytotoxic effects of immune cells on tumors. However, the actual functions of these genes may vary across different cancer types. To characterize the actual functions of sensitizers and resistors within each cancer type, we designed a computational approach (Fig. 3A). Briefly, we first manually curated a list of immune-related features with anti−/pro-tumor activity, including immune cells, immune checkpoint genes, and inflammatory cytokines (Table S5). Then, the associations between the dichotomous status of sensitizers/resistors and the abundance of immune-related features were calculated using a regression-based approach. Sensitizer-related (S-related) features were defined as anti-tumor features with negative association (lower activity upon inactivation) or pro-tumor features with a positive association (higher activity upon inactivation) with sensitizer genes. Resistor-related (R-related) features were defined in an analogous way. Accordingly, the number of S−/R-related features for each sensitizer or resistor genes in each cancer type could be determined (Fig. S5A-B). Conceptually, a sensitizer/resistor with greater S−/R-related feature number is more likely to have a crucial role in anti-tumor immunity. We also calculated the total S−/R-related feature number across different cancers of sensitizers and resistors, respectively (Fig. 3B, C). Notably, it could be observed that some well-recognized immune regulators, such as B2M and TAP1, were relatively top-ranked. Besides, we also found that resistor genes have more related features than sensitizer genes (Fig. 3D).
Identification of MON2 as novel immuno-oncology target
To investigate whether the functional status of sensitizers and resistors was associated with the survival outcome of cancer patients, we conducted Cox proportional hazards regression analysis for all the 105 sensitizer and resistor genes, controlling for confounding factors including cancer type and age. The immune subtype of C2 was found to have an immune-inflamed phenotype, and thus there might exist more interactions between cancer cells and cytotoxic immune cells in tumors of this subtype [25]. Therefore, for the C2 subtype, the impact of the functional status of sensitizers and resistors on prognosis was more likely to be mediated by antitumor immunity-associated mechanisms. Accordingly, only 2591 cases in the C2 subtype were utilized to conduct this analysis, which yielded 3 and 10 significant prognostic sensitizer and resistor genes, respectively (Fig. 4A). Besides, we explored the prognostic associations of sensitizers and resistors in ICB-treated datasets. A total of 10 datasets from eight different studies were collected, which were then integrated into a metadata set (Table S6). The dichotomous functional status of sensitizer and resistor genes in this set was defined similarly to the TCGA Pan-Cancer cohort. Survival analysis in the metadata set identified nine significant prognostic genes, including five sensitizers and four resistors (Fig. 4A). Interestingly, it could be observed that most of the sensitizers (75%) served as unfavorable prognostic factors (hazard ratio > 1 upon inactivation), while most of the resistors (93%) were associated with favorable prognosis (hazard ratio < 1 upon inactivation), which was in accordance with their own properties.
Considering that the loss of a certain gene in itself may affect the survival of cancer cells without the involvement of antitumor immune response, we performed an additional filtering procedure using gene dependency data from the Cancer Dependency Map (DepMap) portal (depmap.org) to avoid potential bias introduced by essential genes [27]. The average CERES scores across 739 cancer cell lines were calculated, and genes with CERES scores ranging from − 0.25 to 0.25 were considered tumor proliferation-independent genes (n = 14,182) (Fig. 4B). By intersecting the above results, a resistor named MON2 was identified as the only proliferation-independent gene, which had significant prognostic associations both in the TCGA C2 and combined ICB-treated cohorts (Fig. 4C).
Association of MON2 with antitumor immunity
Systematic analyses were conducted to characterize the roles of MON2 in antitumor immunity. MON2 was identified as a significant resistor gene by two of the 17 screens (Fig. S6A). Given that one of the screens was carried out in cells derived from mouse breast cancer, we validated our finding in two human breast tumor cell lines, MDA-MB-231 and MCF7. MON2-deficient tumor cells were generated by CRISPR/Cas9. The knockout efficiency was determined by the Tracking of Indels by Decomposition (TIDE) [28] (Fig. S6B). The resulting cells were loaded with a defined antigen and co-cultured with TCR-transduced T cells that recognize the antigen-MHC class I complex. We observed that both MDA-MB-231 and MCF7 cells with MON2 deficiency were more sensitive to the killing effect of T cells than MON2-proficient cells, consistent with the results from the CRISPR screen in mouse cells (Fig. 4D, E).
We further interrogated the relationship between MON2 and the three above scores, including immune, CYT, and MHC scores. MON2 was most significantly associated with MHC scores (Fig. S7A). Since MHC scores were used to measure the antigen presentation ability of tumor cells, we hypothesized that MON2 might be more relevant to tumor cell-intrinsic immune-associated factors. To validate this hypothesis, the associations between MON2 and the other three tumor cell-intrinsic factors, including mutational load, single nucleotide variant (SNV)-based neoantigen load, and insertion and deletion (indel)-based neoantigen load, were tested. As expected, a high association was found between MON2 and all three factors (Fig. S7B). These analyses collectively showed that the loss of MON2 was closely related to enhanced tumor immunogenicity.
Based on the response data from ICB-treated datasets, we next set out to evaluate whether a direct correlation was present between MON2 and response to ICBs. Unsurprisingly, the loss of MON2 was associated with a more favorable response to ICBs in multiple datasets, albeit not very remarkable (Fig. 4F). As a complementary investigation, the relationship between the expression of MON2 and the objective response rate (ORR) of anti-PD-1/PD-L1 therapy was investigated as well [29]. A marginally significant negative correlation between MON2 expression and ORR was obtained (P = 0.093), which could corroborate the above conclusion to some extent (Fig. 4G). Given that MON2 expression could be a potential factor affecting immunotherapy response, we delineated the distribution of MON2 expression across different cancer types (Fig. S7C) and immune subtypes (Fig. S7D). Notably, it could be observed that MON2 exhibited the highest expression level in the C5 subtype (immunologically quiet). Furthermore, the crucial role of microsatellite instability (MSI) status in the field of cancer immunotherapy prompted us to examine the association between MON2 expression and MSI status [30]. Three cancer types in TCGA, including COAD, STAD, and UCEC, were selected for conducting this analysis with a large sample size of MSI-high (MSI-H) tumors [31]. MSI-low (MSI-L)/microsatellite-stable (MSS) tumors showed a significantly higher expression of MON2 in COAD, STAD, and the combined cohorts (Fig. S8A, B, D), except for UCEC (Fig. S8C). Correlation analysis between MSI events and MON2 expression was also conducted, and a significant negative correlation was obtained, as expected (P = 0.004) (Fig. S8E). Of note, MSI-L/MSS tumors as well as tumors in C5 subtypes were considered to be resistant to currently approved ICBs [25, 32, 33]. Collectively, inhibition of MON2 might lead to the augmentation of anti-tumor immunity, which might open new possibilities for treating tumors with high MON2 activity.
Construction of CTIS for predicting immunotherapy response
Despite recent advances, predicting response to immunotherapy remains challenging and requires further investigations. As our findings suggested important roles of sensitizer and resistor genes in immune cell-mediated tumor killing, we postulated that it might be clinically meaningful to construct a predictive signature based on these genes.
The detailed process of signature construction was displayed in Fig. 5A. Only pretreatment samples were included in this step. A preliminary filtering was first performed to exclude genes with less than 100 S−/R-related features (across all the cancer types). As a result, 39 sensitizers and 37 resistors were retained after the filtration. Then, based on Liu et al.’s dataset which has the largest pretreatment sample size, minimal depth (MD)-based random survival forest (RSF) analysis was conducted to further narrow down the gene list [34]. The RSF analysis was repeated 1000 times and six genes that led to a largest concordance index (C-index) value was considered as the final candidates. These genes consisted of four sensitizers (JAK1, NFKB2, PPP6C and TNFRSF1B) and two resistors (PIGM and TPR). A CRISPR screening-based tumor-intrinsic immune score (CTIS) was then defined as the difference between the mean expression of four sensitizers and two resistors and a higher CTIS indicated a stronger anti-tumor immune response. In the discovery dataset (Liu et al.’s dataset), CTIS exhibited a significant prognostic value (P < 0.0001); higher CTIS was related to a better survival outcome (Fig. 5B). We also investigated its performance in several validation datasets. Among them, we found CTIS was significantly associated with prognosis in two datasets and patients with higher CTIS consistently had an improved overall survival (Fig. 5C).
After that, we comprehensively compared the CTIS signature with other 14 published signatures (details in Supplemental Materials and Methods). The association between these signatures was first investigated using the expression data of ICB-treated datasets. It could be observed that many signatures, such as CD8, CYT, IFNG, PD-L1, and inflamed signatures, were highly correlated with each other (Fig. 5D). In contrast, the CTIS signature presented relatively moderate associations with other signatures, suggesting it has a complementary rather than an alternative role. Then, comparison of the performance for predicting response to ICB treatment between CTIS and other signatures was performed (Fig. S9). Although CTIS was not the best biomarker across all the datasets, it outperformed other signatures in the condition of comparing the overall performance, with a mean AUC value of 0.620 (Fig. 5E).
We further calculated the CTIS of treatment-naïve tumors from TCGA Pan-Cancer cohort to characterize the distribution of CTIS across different cancer types (Fig. S10A) and immune subtypes (Fig. S10B). Among the six immune subtypes, C5 had the lowest CTIS scores, consistent with its property. Although higher CTIS indicated a better prognosis in ICB-treated datasets, in treatment-naïve TCGA tumors, CTIS conferred a dual prognostic impact depending on the cancer type (Fig. S10C). In addition, relationships among mutation and neoantigen load and CTIS were also assessed. However, no remarkable results were observed (Fig. S10D).
Determination of potential OGs/TSGs regulatory network
It has been widely reported that some oncogenes (OGs) or tumor suppressor genes (TSGs) may be involved in the regulation of tumor immunity [35]. For example, activation of MYC can induce the transcription of both CD47 and PD-L1 and thereby suppress the anti-tumor immune response [36]. Accordingly, we considered that it was reasonable to build an OGs/TSGs regulatory network for sensitizers and resistors, which could uncover potential mechanisms underlying the dysregulation of sensitizers and resistors. The list of OGs and TSGs was obtained from OncoKB (oncokb.org) [37]. For OGs, functionally relevant events were considered as activation (0 = non-activation; 1 = activation). As for TSGs, functionally relevant events were considered as inactivation (0 = non-inactivation; 1 = inactivation). We hypothesized that associations existed between sensitizers and OGs/TSGs when the expression of sensitizers was significantly down-regulated upon OG activation/TSG inactivation (adj. P < 0.05 and log-fold change < − 0.25). Correspondingly, an association between resistors and OGs/TSGs was found when resistors’ expression was significantly up-regulated upon OG activation/TSG inactivation (adj. P < 0.05 and log-fold change > 0.25). We identified a total of 159 significant associations between sensitizers/resistors and OGs/TSGs (Fig. 6A). Interestingly, we found that sensitizer genes tended to have more associations with TSGs (66.7%), while resistor genes exhibited considerably more associations with OGs (86.4%) (Fig. 6B). OGs/TSGs with significant associations were defined as tumor immunity-related OGs/TSGs (Table S7).
Prediction of repurposing candidates synergistic with immunotherapies
One major challenge in cancer immunotherapy is to increase tumor response to ICB treatment. Previous experimental and clinical studies have demonstrated that combinatorial therapeutic strategies could substantially increase the percentage of responder cases and contribute to significant survival benefits [38, 39]. A signature matching approach for drug prediction was adopted to discern more potential combination partners synergistic with immunotherapies (Fig. 6C). Since query signature was the basis of this approach, we first sought to collect potential immunotherapy-relevant genes for drug retrieval. Data from additional five ICB-treated screens were obtained. The experimental design used in these screens differed from ICB-naïve screens, which also introduced immunotherapeutic intervention into the experiments (Fig. S11A). Therefore, these screens could identify potential regulators mediating resistance (ICB enhancer genes) or sensitivity (ICB suppressor genes) to ICBs upon loss (Table S8). These screens studied four cancers (Fig. S11B) and common hits of enhancers and suppressors were also displayed, respectively (Fig. S11C). Since similar properties existed between sensitizers and enhancers as well as between resistors and suppressors, these genes with intersected (Fig. S11D, E). There were multiple common genes between the lists of sensitizers and enhancers (Fig. S11D). Subsequently, sensitizer genes, tumor immunity-related TSGs, and ICB enhancer genes were integrated into a meta-gene list, and genes in this list were termed positive regulators. Theoretically, potentiating the function of these genes might be related to an enhanced anti-tumor immune response. The list of negative regulators was obtained in the same way. Positive and negative regulators together constituted the query signature (Fig. 6c, Table S9).
The other necessary component of the signature matching-based approach was the drug signatures, also known as drug-induced profiles of expression changes. These drug signatures were downloaded from the Connectivity Map (CMap) datasets (CMap Build 2: 1288 compounds) [40]. To ensure easy clinical translation of our findings, we only selected drugs that have been approved or have already passed phase I and II clinical trials for subsequent analysis, leveraging the annotations from the Drug Repurposing Hub [41]. To match query signature with drug signatures, we applied three methods, including eXtreme Sum (XSum) [42], Kolmogorov-Smirnov (KS) [40], and the Reverse Gene Expression Score (RGES) [43]. Due to the varying measurement scales of the resultant scores from the three methods, an order statistics-based method developed by Stuart et al. was utilized to integrate these results and yielded a robust drug prediction result (Fig. 6C) [44]. The results of drug prediction are presented in Fig. 6D and Table S10. Notably, among the top 10 drugs, four candidates, including irinotecan, quercetin, trifluridine, and resveratrol, were previously reported by experimental or clinical studies, demonstrating the reliability of our prediction.
