Research Paper Volume 15, Issue 7 pp 2667—2688

Machine learning-based construction of immunogenic cell death-related score for improving prognosis and response to immunotherapy in melanoma

Guoyin Li1,2,3, , Huina Zhang1, , Jin Zhao1, , Qiongwen Liu1, , Jinke Jiao1, , Mingsheng Yang1, , Changjing Wu1, ,

  • 1 College of Life Science and Agronomy, Zhoukou Normal University, Zhoukou, Henan, China
  • 2 Key Laboratory of Modern Teaching Technology, Ministry of Education, Shaanxi Normal University, Xi’an, Shaanxi, China
  • 3 Academy of Medical Science, Zhengzhou University, Zhengzhou, Henan, China

Received: January 8, 2023       Accepted: March 24, 2023       Published: April 6, 2023      

https://doi.org/10.18632/aging.204636
How to Cite

Copyright: © 2023 Li et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

Background: Immunogenic cell death (ICD) is a form of regulated cell death (RCD) which could drive the activation of the innate and adaptive immune responses. In this work, we aimed to develop an ICD-related signature to facilitate the assessment of prognosis and immunotherapy response for melanoma patients.

Methods: A set of machine learning methods, including consensus clustering, non-negative matrix factorization (NMF) method and least absolute shrinkage and selection operator (LASSO) logistic regression model, and bioinformatics analytic tools were integrated to construct an ICD-related risk score (ICDscore). CIBERSORT and ESTIMATE algorithm were used to evaluate the infiltration of immune cells. The 'pRRophetic' package in R and 6 cohorts of melanoma patients receiving immunotherapy were used for therapy sensitivity analyses. The predictive performance between ICDscore with other mRNA signatures were also compared.

Results: The ICDscore could predict prognosis and immunotherapy response in multiple cohorts, and displayed superior performance than other forms of cell death-related signatures or 52 published signatures. The melanoma patients with low ICDscore were marked with high infiltration of immune cells, high expression of immune checkpoint inhibitor-related genes, and increased tumor mutation burden.

Conclusions: In conclusion, we constructed a stable and robust ICD-related signature for evaluating the prognosis and benefits of immunotherapy, and it could serve as a promising tool to guide decision-making and surveillance for individual melanoma patients.

Introduction

An ideal anti-cancer therapy is still an ongoing pursuit for most researchers and clinicians. The extensive research efforts in past decades have resulted in the emergence of immune-checkpoint inhibitors (ICIs) and various types of targeted therapy. However, not all patients, even with the same tumor type and at the same stage, responded well to these anti-tumor agents. For instance, only 4% uveal melanoma patients with metastasis showed partial response to nivolumab plus ipilimumab, as reported in a recent published real-life, retrospective study [1]. In a multicenter, randomized phase 3 study, only 36-37% and 13% of melanoma patients at advanced stage showed an objective response to pembrolizumab and ipilimumab, respectively [2]. Given the complexity of the human genome and the heterogeneity of the tumor itself, personalized therapy might be an alternative option to deal with the aforementioned problem. Some efforts have been made to classify patients into subgroups with a different risk score, based on the generation of gene expression signatures [3]. The most successful signature might be the 21-gene expression assay, which helps to identify a subgroup of hormone-receptor–positive, HER2-negative, axillary node–negative breast cancer patients with a high risk of recurrence and to guide adjuvant chemotherapy for these patients [4, 5].

Inducing the death of tumor cells is the ultimate purpose of cancer treatment. Multiple forms of regulated cell death (RCD) have been identified in cancer, including apoptosis, autophagy, ferroptosis, necroptosis, and proptosis [6]. Immunogenic cell death (ICD) is one kind of RCD but unique in its capability to trigger antigen-specific adaptive immunological responses [7]. Various types of anti-cancer therapy, including chemotherapy and radiotherapy, could drive the occurrence of ICD [7]. In addition, induction of ICD also becomes a strategy in the design of anti-cancer agents in various types of tumors including melanoma. For example, Li et al. generated polymer micelles (named PPIR780-ZMS) containing IR780 dye and manganese zinc sulfide nanoparticles (ZMS), and revealed that PPIR780-ZMS could maximize ICD in melanoma and augment the infiltration of cytotoxic T cells (CD8+) and helper T cells (CD4+) [8]. Zhang et al. synthesized a nano-inducer (named DOX/ADS NP) via the self-assembly of doxorubicin (DOX) and dermatan sulphate derivative (ADS) grafted by aromatic thioketal (ATK) [9]. This nano-inducer could enhance ICD and activate the immune response in melanoma, leading to shrinkage of tumor [9].

Although some previous studies have attempted to develop signatures, based on some form of RCD, like autophagy and ferroptosis [1013], and classify melanoma patients into subgroups with a different risk score, little is known about the application of ICD-based signature in melanoma. Since melanoma is supposed to be an immunogenic tumor whose proliferation and progression have a tight relationship with immune cells [14, 15], and ICIs have been shown by clinical trials to improve the survival benefits of melanoma patients and have been approved for clinical use globally [1618], we hypothesize that a novel ICD-based risk score might be a better signature in predicting prognosis and immunotherapy efficiency of melanoma patients.

In this work, unsupervised clustering of melanoma patients, based on 34 ICD-related genes summarized by a previous study [19], indicated that these patients could be classified into two subgroups with distinct overall survival (OS). Subsequently, we constructed an ICD-related risk score (ICDscore) and found that melanoma patients with high ICDscore had apparently shorter OS and poor sensitivity to immunotherapies. More importantly, we compared ICDscore with other signatures developed on other forms of RCD and noticed that ICDscore had better performance in predicting OS and immunotherapy efficiency. We also compared ICDscore with 52 previously published signatures, and the results indicated that ICDscore had certain advantages over these signatures regarding prognosis predictability.

Materials and Methods

Public data acquisition and processing

The TCGA-SKCM cohort was downloaded from the Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov). The GSE35640, GSE54467, GSE22153, and GSE65904 were obtained from the Gene Expression Omnibus (GEO) database (https://ncbi.nlm.nih.gov/gds). The gene expression data and clinical information of the Peking University Cancer Hospital (PUCH) study, Riaz17 study, Liu19 study, VanAllen15 study and Gide19 study were downloaded from the GitHub website (https://github.com/), as reported in Cui’s study [20]. All datasets were processed as described in our previous work [21]. All datasets used in this work were downloaded from public databases, an extra ethical approval was not necessary.

Construction of ICDscore

The flow to generate the ICDscore is summarized in Figure 1. Specifically, melanoma patients were firstly stratified into two clusters (namely cluster A and B) via the consensus clustering and non-negative matrix factorization (NMF) clustering methods, based on the gene expression of 34 ICD-related genes. For the consensus clustering, the ‘ConcensusClusterPlus’ package in R was used [22], and the parameters used in the analyses were maxK = 10, reps = 1, 000, pItem = 0.8, pFeature = 1, clusterAlg = “pam”, corUse = “complete.obs”, seed = 123456. For the NMF clustering, the “NMF” package in R was used [23], and the ranks were set from 2 to 10 to do the NMF rank survey. The differentially expressed genes (DEGs) between cluster A and B were assessed by the “limma” package in R. Univariate Cox regression analyses were then performed to identify DEGs which had a significant p-value < 0.1 in the TCGA-SKCM and GSE65904 cohorts. Subsequently, a total of 182 DEGs were then input into a Least absolute shrinkage and selection operator (LASSO) regression model in GSE65904. 4 key genes were generated, and their corresponding coefficients were obtained via multi-variate Cox analysis. The risk score for each patient was calculated by the following formula: score = - 0.17244 * GBP2 - 0.14202 * LYZ - 0.24610 * CST7 - 0.19525 * SIRPG. The ICDscore of patients in each cohort was calculated with the formula: ICDscore = (score-Min) / absolute (Max), as reported in our previous studies [21, 24].

The work flow for the construction of ICDscore. Two clustering methods (consensus clustering and NMF clustering) were used for the molecular subtyping of melanoma patients, based on the gene expression of 34 ICD related genes. Two clusters (named as Cluster A and B) were identified and DEGs between these two clusters were analyzed. The LASSO regression model, multivariate Cox analyses were then used for the construction of ICDscore. The association between the ICDscore and prognosis, tumor immune microenvironment, or immunotherapy response was comprehensively investigated. The performance between ICDscore and other signatures was compared.

Figure 1. The work flow for the construction of ICDscore. Two clustering methods (consensus clustering and NMF clustering) were used for the molecular subtyping of melanoma patients, based on the gene expression of 34 ICD related genes. Two clusters (named as Cluster A and B) were identified and DEGs between these two clusters were analyzed. The LASSO regression model, multivariate Cox analyses were then used for the construction of ICDscore. The association between the ICDscore and prognosis, tumor immune microenvironment, or immunotherapy response was comprehensively investigated. The performance between ICDscore and other signatures was compared.

Immune profile analysis

The infiltration ratio of 22 immune cells in patients was calculated by the CIBERSORT algorithm in R software, as reported in our previous study [21]. The immunescore and stromalscore of each patient were calculated by the ‘estimate’ package in R [25].

Enrichment analysis

Gene Set Enrichment Analysis (GSEA) of SKCM patients were performed by the ‘clusterProfiler’ package in R. The c5.go.bp.v2022.1.Hs.symbols.gmt was chosen as the gene set database. The ‘GseaVis’ package in R was used for visualization [26].

Statistical analysis

All the data were processed, analyzed and visualized by R software (version 4.1.3). In addition to the packages mentioned above, other packages in R used in this work included “tidyverse”, “survival”, “msigdbr”, “dplyr”, “org.Hs.eg.db”, “ggplot2”, “glmnet”, “scales”, “aplot ”, “survivalROC”, “ggrepel”, “enrichplot”, “corrplot”, “survminer”, “timeROC”, “rms”, “pec”, “ggalluvial”, “VennDiagram”, “ggh4x”, “patchwork”, “pRRophetic”, and “CompareC”. The Kaplan-Meier method was used for prognosis analyses. The Correlation analyses were conducted with the Pearson method. The comparison of categorical variables between two groups was conducted with the chi-square test. The continuous variables were compared with the Wilcoxon rank-sum test. A value of p < 0.05 was considered to be statistically significant (*, p < 0.05; **, p < 0.01; ***, p < 0.001; ****, p < 0.0001).

Results

Unsupervised clustering of ICD-related genes in melanoma

Based on the expression of 34 ICD-related genes identified by Abhishek D. Garg et al. [19], two unsupervised clustering methods were employed to stratify melanoma patients in the TCGA-SKCM cohort. For the Consensus Clustering method, k = 2 was selected as the optimal parameter for further analyses, based on the consensus matrix for k = 2 to 10 and the consensus cumulative distribution function (CDF) plot (Figure 2A, 2B and Supplementary Figure 1A1H). Two hundred thirty-five melanoma patients were classified into C1 cluster and showed significantly prolonged median overall survival (OS) than those in the C2 cluster (Figure 1C). For the Nonnegative Matrix Factorization (NMF) method, rank = 2 was chosen as the most suitable number of subgroups according to the cophenetic coefficient and the consensus matrix for different rank numbers (Figure 2D, 2E and Supplementary Figure 1I, 1J). 194 patients were stratified into C1 cluster and exhibited better prognosis (Figure 2F, p < 0.0001). No matter which method was applied, an obviously different distribution between was observed the C1 and C2 clusters via the PCA plots (Figure 2G and Supplementary Figure 1K). As shown in Figure 2H, 192 melanoma patients were stratified in the C1 cluster by both methods and named Cluster A. Meanwhile, 210 patients were grouped in the C2 cluster by both methods and were named as Cluster B. The median OS of melanoma patients in Cluster A reached 4634 days, and was considerably longer than that of patients in Cluster B (1766 days, Supplementary Figure 1L, p < 0.0001). In addition, the transcriptional expression of most of the ICD-related genes was significantly elevated in patients of Cluster A (Figure 2H).

Clustering of melanoma patients based on ICD related genes. (A) The cumulative distribution function (CDF) curve of consensus clustering for k = 2 to 10. (B) Heatmap depicts consensus clustering solution (k = 2) for 34 ICD related genes in melanoma patients from the TCGA-SKCM dataset. (C) Kaplan–Meier curves of OS in the C1 and C2 subtypes (the consensus clustering method) of melanoma patients. (D) The optimal rank was selected as 2 since the cophenetic coefficient firstly started decreasing at this point. (E) Heatmap of NMF clustering results of melanoma patients from the TCGA-SKCM dataset. (F) Kaplan–Meier curves of OS in the C1 and C2 subtypes (the NMF method) of melanoma patients. (G) Principal component analysis (PCA) on the expression level of 34 ICD related genes in clusters classified by NMF method. (H) Venn diagram to identify melanoma patients in the C1 cluster (Cluster A) and C2 cluster (Cluster B) defined by both clustering methods. (I) Gene expression comparison of 34 ICD related genes between Cluster A and B in the TCGA-SKCM cohort. Ns, not significant; *p

Figure 2. Clustering of melanoma patients based on ICD related genes. (A) The cumulative distribution function (CDF) curve of consensus clustering for k = 2 to 10. (B) Heatmap depicts consensus clustering solution (k = 2) for 34 ICD related genes in melanoma patients from the TCGA-SKCM dataset. (C) Kaplan–Meier curves of OS in the C1 and C2 subtypes (the consensus clustering method) of melanoma patients. (D) The optimal rank was selected as 2 since the cophenetic coefficient firstly started decreasing at this point. (E) Heatmap of NMF clustering results of melanoma patients from the TCGA-SKCM dataset. (F) Kaplan–Meier curves of OS in the C1 and C2 subtypes (the NMF method) of melanoma patients. (G) Principal component analysis (PCA) on the expression level of 34 ICD related genes in clusters classified by NMF method. (H) Venn diagram to identify melanoma patients in the C1 cluster (Cluster A) and C2 cluster (Cluster B) defined by both clustering methods. (I) Gene expression comparison of 34 ICD related genes between Cluster A and B in the TCGA-SKCM cohort. Ns, not significant; *p < 0.05; **, p < 0.01; ***, p < 0.001; ****, p < 0.0001.

The biological features of melanoma patients in the two clusters

Besides ICD, various other forms of cell death have been identified, such as autophagy, ferroptosis and necroptosis [27]. These forms of cell death have some shared regulators or signaling components, but also have unique biological characteristics. Some previous studies tried to classify melanoma patients into high- and low-risk subgroups based on the features of these forms of cell death [1013]. As shown in Figure 3A, most the melanoma patients in Cluster A were labeled as low risk when they were classified based on the ferroptosis-, autophagy-, proptosis-, or necroptosis-related signatures. The concordance between ICD-based classification and other forms of cell death-based classification was 75.62% (ferroptosis), 58.71% (autophagy), 74.38% (proptosis), and necroptosis (84.33%). To characterize the biological features of patients in Cluster A and Cluster B, GSEA was employed. The results showed that melanoma patients in Cluster A showed enrichment in immune related pathways (Supplementary Table 1) like lymphocyte mediated immunity (Figure 3B) and activation of immune response (Figure 3C). On the other hand, patients in the Cluster B showed an enrichment in the processing and translation of RNA (Supplementary Table 1 and Figure 3D, 3E), and in DNA repair (Supplementary Table 1 and Figure 3F). Consistently, melanoma patients had an obviously higher level of infiltration of anti-tumor immune components, such as CD8 T cells, M1 macrophage, and activated memory CD4 T cells (Figure 3G), whereas some pro-tumor immune components, like M2 macrophage, resting mast cells, and resting memory CD4 T cells, showed a significantly higher infiltration in patients from the Cluster B (Figure 3G).

Biological features of melanoma patients in Cluster A and B. (A) Sankey diagram showed the connection degree between ICD-based classification and other forms of cell death-based classification. (B–F) Examples of GSEA results of melanoma patients in Cluster A (B, C) and B (D–F). (G) Distribution of 22 types of infiltrating immune cells in melanoma patients in the Cluster A and B. Ns, not significant; *p

Figure 3. Biological features of melanoma patients in Cluster A and B. (A) Sankey diagram showed the connection degree between ICD-based classification and other forms of cell death-based classification. (BF) Examples of GSEA results of melanoma patients in Cluster A (B, C) and B (DF). (G) Distribution of 22 types of infiltrating immune cells in melanoma patients in the Cluster A and B. Ns, not significant; *p < 0.05; **, p < 0.01; ***, p < 0.001; ****, p < 0.0001.

Construction of ICD-related signature in melanoma

Unsupervised clustering revealed that ICD-related genes could help to stratify melanoma patients into subgroups with distinct prognoses and tumor microenvironment (TME) (Figures 2, 3). To construct a signature that help to recognize these two subgroups, we first analyzed the differentially expressed genes (DEGs) between the patients in Cluster A and those in Cluster B. A total of 1, 038 DEGs were identified in the TCGA-SKCM cohort with the criteria of absolute logFC (fold change) ≥ 1 and adjusted p-value < 0.05 (Supplementary Table 2). To better reflect the difference between these two clusters, only DEGs (genes marked with red, Supplementary Table 2) with absolute logFC ≥ 2 were selected for further analyses. The subsequent univariate Cox analyses in TCGA-SKCM and GSE65904 indicated that 182 DEGs had a significant p-value < 0.1 in both datasets (Figure 4A). The 182 DEGs were then input into a LASSO regression model in GSE65904 as described in the Method section (Figure 4B, 4C). Four crucial genes were generated, and they were guanylate binding protein 2 (GBP2), lysozyme (LYZ), cystatin F (CST7), and signal regulatory protein gamma (SIRPG). The ICDscore was calculated based on the transcriptional expression and coefficient of these 4 genes, as demonstrated in the Method section. When melanoma patients were divided into two groups based on the median value of ICDscore in each cohort, those patients in the high-ICDscore subgroup showed a significantly shorter OS in the training (GSE65904, Figure 4D) and validating (TCGA-SKCM, GSE54467, and GSE22153, Figure 4E4G) datasets. The forest plot of meta-analysis (Supplementary Figure 2A) also indicated ICDscore as a vital risk factor for melanoma patients. Besides, melanoma patients with high ICDscore also exhibited considerably shorter progression-free survival (PFS, Supplementary Figure 2B, 2C). As shown in Figure 4H, patients with low ICDscore had apparently higher expression of all the four genes. In addition, the high-ICDscore subgroup had a significantly higher percentage of melanoma patients with deeper Breslow depth (p < 0.01), advanced Clark level (p < 0.001), advanced T stage (p < 0.01), and dead status (p < 0.0001). Correspondingly, melanoma patients with deeper Breslow depth (Supplementary Figure 2G), advanced Clark level (Supplementary Figure 2H), advanced T stage (Supplementary Figure 2I) exhibited a significantly higher level of ICDscore, and no difference of ICDscore level was observed in melanoma patients divided by gender (Supplementary Figure 2E), age (Supplementary Figure 2F), N stage (Supplementary Figure 2J), M stage (Supplementary Figure 2K) or clinical stage (Supplementary Figure 2L). Besides, univariate and subsequent multivariate Cox analyses revealed that ICDscore could serve as an independent prognostic factor (Figure 4I, 4J).

Construction of ICDscore. (A) Venn diagram to screen DEGs with a significant prognostic p-value B, C) The LASSO Cox regression model was constructed from 182 DEGs, and the tuning parameter (λ) was calculated based on the partial likelihood deviance with ten-fold cross validation. 4 signature genes were identified according to the best fit profile. (D–G) Kaplan–Meier curves of OS in melanoma patients from ICDscore-high and ICDscore-low subgroups of GSE65904 (D), TCGA-SKCM (E), GSE54467 (F), and GSE22153 (G) datasets. (H) Clinical characteristics and RNA expression level of 4 crucial genes in melanoma patients from ICDscore-high and ICDscore-low subgroups of the TCGA-SKCM dataset. (I) Univariate analysis of ICDscore and other clinical characteristics in TCGA-SKCM dataset. (J) Multivariate analysis shows ICDscore, breslow depth, M stage and age were independent prognostic factors. Ns, not significant; *p

Figure 4. Construction of ICDscore. (A) Venn diagram to screen DEGs with a significant prognostic p-value < 0.1 in the TCGA-SKCM and GSE65904 cohorts. (B, C) The LASSO Cox regression model was constructed from 182 DEGs, and the tuning parameter (λ) was calculated based on the partial likelihood deviance with ten-fold cross validation. 4 signature genes were identified according to the best fit profile. (DG) Kaplan–Meier curves of OS in melanoma patients from ICDscore-high and ICDscore-low subgroups of GSE65904 (D), TCGA-SKCM (E), GSE54467 (F), and GSE22153 (G) datasets. (H) Clinical characteristics and RNA expression level of 4 crucial genes in melanoma patients from ICDscore-high and ICDscore-low subgroups of the TCGA-SKCM dataset. (I) Univariate analysis of ICDscore and other clinical characteristics in TCGA-SKCM dataset. (J) Multivariate analysis shows ICDscore, breslow depth, M stage and age were independent prognostic factors. Ns, not significant; *p < 0.05; **, p < 0.01; ***, p < 0.001; ****, p < 0.0001.

Immune landscape of ICDscore-stratified melanoma patients

As shown in Supplementary Figure 2C, 94.79% of patients in the Cluster A were classified into low-ICDscore subgroup, while 90.52% of patients in the Cluster B were into the high-ICDscore subgroup. Most of ICD-related genes showed a significantly elevated expression in patients from the low-ICDscore group (Figure 5A). Correlation analyses in both the TCGA-SKCM and GSE65904 cohorts revealed that ICDscore had a strong negative correlation with CD8 T cells (R = -0.65 in TCGA-SKCM cohort and -0.53 in GSE65904 cohort, Figure 5B), and significantly positive correlation with NK cells resting, M0 macrophage, M2 macrophage, or mast cell resting (Figure 5B). In addition, melanoma patients from the ICDscore-low subgroup had an apparently higher immunescore and stromalscore (Figure 5C, 5D), suggesting a higher infiltrative level of immune and stromal cells.

Immune profile of ICDscore-based classification. (A) Gene expression comparison of 34 ICD related genes between ICDscore-high and ICDscore-low subgroups in the TCGA-SKCM cohort. (B) Correlation analyses between ICDscore and infiltration level of 22 immune cells in the TCGA-SKCM and GSE65904 datasets. (C, D) Comparison of immunescore (C) and stromalscore (D) between ICDscore-high and ICDscore-low subgroups in the TCGA-SKCM and GSE65904 datasets. (E) Box plot showing a difference in the value of ICDscore across the five subtypes for melanoma patients in the TCGA

Figure 5. Immune profile of ICDscore-based classification. (A) Gene expression comparison of 34 ICD related genes between ICDscore-high and ICDscore-low subgroups in the TCGA-SKCM cohort. (B) Correlation analyses between ICDscore and infiltration level of 22 immune cells in the TCGA-SKCM and GSE65904 datasets. (C, D) Comparison of immunescore (C) and stromalscore (D) between ICDscore-high and ICDscore-low subgroups in the TCGA-SKCM and GSE65904 datasets. (E) Box plot showing a difference in the value of ICDscore across the five subtypes for melanoma patients in the TCGA_SKCM dataset. (F) Sankey diagram showed the connection degree between ICDscore-based classification and Akbani cluster, TME subtype and mutation subtype in the TCGA-SKCM dataset. (G) Box plot showing a difference in the expression of multiple exhausted T cell markers or immune-checkpoint markers between ICDscore-high and ICDscore-low subgroups in the TCGA-SKCM cohort. (H) Box plot showing a difference in the TMB between ICDscore-high and ICDscore-low subgroups in the TCGA-SKCM cohort. Ns, not significant; *p < 0.05; **, p < 0.01; ***, p < 0.001; ****, p < 0.0001.

The relationship between ICDscore-based classification and other method-based classification was also investigated. TCGA research network divided cancer patients into six clusters [28]. As shown in Figure 5E, the C2 cluster (INF-gamma dominant) had the lowest ICDscore, whereas the C4 (lymphocyte depleted) and C1 (wound healing) clusters had relatively high ICDscore. Besides, most melanoma patients in the low-ICDscore subgroup were categorized into the immune cluster, which was defined on the basis of consensus hierarchical clustering of 1500 genes [29]. Similarly, 83.26% of patients in the low-ICDscore subgroup were labeled as ‘Immune-Enriched, Fibrotic’ (IE/F) or ‘Immune-Enriched, Non-Fibrotic’ (IE), based on the characterization of TME by 29 functional gene expression signatures (Fges) [30]. In addition, melanoma patients were also divided into four subtypes based on some mutated genes, and they were B-Raf Proto-Oncogene, Serine/Threonine Kinase (BRAF) hotspot, RAS hotspot, Neurofibromin 1 (NF1) mutant, and Triple-WT (wild-type) [29]. The distribution of ICDscore-high or ICDscore-low patients in BRAF hotspot, NF1 mutant, or RAS hotspot subtype was similar (Figure 5F, Chi-square test, p > 0.05), but a significantly higher ratio of ICDscore-high patients was Triple-WT (Figure 5F, Chi-square test, p = 0.0465).

Meanwhile, the expression of multiple exhausted T cell markers or immune-checkpoint markers, including CD27, CD274 (also referred as PDL1), cytotoxic T-lymphocyte associated protein 4 (CTLA4), programmed cell death 1 (PDCD1), hepatitis A virus cellular receptor 2 (HAVCR2, also named as TIM3), T cell immunoreceptor with Ig and ITIM domains (TIGIT), thymocyte selection associated high mobility group box (TOX), TNF receptor superfamily member 9 (TNFRSF9), and ectonucleoside triphosphate diphosphohydrolase 1 (ENTPD1) [31, 32], were dramatically elevated in patients with low ICDscore (Figure 5G), suggesting these patients might be able to benefit from immune checkpoint inhibitors (ICIs). Further, tumor mutation burden (TMB), another biomarker in predicting immunotherapy efficiency [33], was also significantly higher in the ICDscore-low subtype (Figure 5H, p < 0.05).

Evaluation of anti-tumor therapy in ICDscore-based subgroups

To support the hypothesis that ICDscore might help to predict immunotherapy efficiency, several cohorts of melanoma patients receiving various forms of immunotherapy were investigated. As shown in Figure 6A, melanoma patients who responded to adoptive T-cell therapy (ACT) treatment (GSE35640), anti-PD-1 (PUCH, Gide19, and Liu19 cohorts) or anti-CTLA-4 (VanAllen15 cohort) showed a significant lower ICDscore than those who were non-responders to immunotherapies. In the Riaz17 cohort, melanoma patients who responded to nivolumab (anti-PD-1 agent) had lower ICDscore than the non-responders, but the difference did not reach significance (Figure 6A). When melanoma patients were divided into ICDscore-high and ICDscore-low subgroups by the median value of ICDscore in each cohort, the patients in the ICDscore-low subgroup showed a higher object responsive rate (ORR) to immunotherapies, and the difference reached significance in the PUCH (p = 0.037), Gide19 (p = 0.018), and Liu19 (p = 0.014) cohorts (Figure 6B). In four cohorts with survival information available, patients in the ICDscore-low subgroup showed apparently longer OS than those in the ICDscore-high subgroup (Figure 6C6F).

Sensitivity evaluation to anti-cancer therapy. (A) ICDscore of melanoma patients receiving immunotherapy in GSE35640, PUCH, Gide19, Liu19, VanAllen16, and Riaz17 cohorts. (B) Ratio of patients responding or not responding to immunotherapy in ICDscore-high and ICDscore-low subgroups of GSE35640, PUCH, Gide19, Liu19, VanAllen16, and Riaz17 cohorts. (C–F) Kaplan–Meier curves of OS in melanoma patients from ICDscore-high and ICDscore-low subgroups of VanAllen15 (C), PUCH (D), Liu19 (E), and Gide19 (G) cohorts. (G–N) Box plot showing a difference in the IC50 values of STF-62247 (G), Erlotinib (H), Rapamycin (I), Sunitinib (J), CGP-60474 (K), JW-7-52-1 (L), Roscovitine (M), CAL-101 (N) between ICDscore-high and ICDscore-low melanoma patients. Ns, not significant; *p

Figure 6. Sensitivity evaluation to anti-cancer therapy. (A) ICDscore of melanoma patients receiving immunotherapy in GSE35640, PUCH, Gide19, Liu19, VanAllen16, and Riaz17 cohorts. (B) Ratio of patients responding or not responding to immunotherapy in ICDscore-high and ICDscore-low subgroups of GSE35640, PUCH, Gide19, Liu19, VanAllen16, and Riaz17 cohorts. (CF) Kaplan–Meier curves of OS in melanoma patients from ICDscore-high and ICDscore-low subgroups of VanAllen15 (C), PUCH (D), Liu19 (E), and Gide19 (G) cohorts. (GN) Box plot showing a difference in the IC50 values of STF-62247 (G), Erlotinib (H), Rapamycin (I), Sunitinib (J), CGP-60474 (K), JW-7-52-1 (L), Roscovitine (M), CAL-101 (N) between ICDscore-high and ICDscore-low melanoma patients. Ns, not significant; *p < 0.05; **, p < 0.01.

We further manipulated the ‘pRRophetic’ package in R software to estimate the drug sensitivity of melanoma patients in ICDscore-classified subgroups. As shown in Supplementary Table 3, ICDscore had a strong correlation (absolute R > 0.5, p-value < 0.05) with patients’ sensitivity to 34 anti-cancer agents. The top 8 agents, according to the coefficiency, were JW-7-52-1, Roscovitine (CDK inhibitor), Rapamycin (mTOR inhibitor), CGP-60474 (CDK inhibitor), Erlotinib (EGFR inhibitor), CAL-101 (PI3K inhibitor), STF-62247 (autophagy inducer), Sunitinib (Multi-kinase inhibitor), and melanoma patients with high ICDscore had a dramatically higher IC50 value than those with low ICDscore (Figure 6G6N).

Comparison of ICDscore and other gene expression-based prognostic signatures

As shown in Figure 3A, unsupervised clustering of melanoma patients based on ICD-related genes had 58.71% to 84.33% concordance with other forms of cell death-based classification. The C-index [95% confidence interval] of ICDscore was 0.669 [0.641 – 0.697], 0.614 [0.592 – 0.636], 0.690 [0.649 – 0.732], and 0.648 [0.605 – 0.692] in GSE65904, TCGA-SKCM, GSE54467 and GSE22153 cohort, respectively. When compared with clinical features of melanoma patients, the C-index of ICDscore was lower than that of Breslow depth, T stage, and Clark level, but higher than that of the rest clinical features (Supplementary Figure 2M). However, no significance was observed in the comparison of ICDscore and clinical features of melanoma patients (Supplementary Figure 2M). As displayed in Figure 7A7D, ICDscore showed certain improved accuracy than ferroptosis-, autophagy-, proptosis-, or necroptosis-based signature (except for comparison between ICDscore and ferroptosis-related risk score in TCGA-SKCM, Figure 7B). Besides, the AUC of ICDscore in predicting efficiency to immunotherapies was 0.701, 0.689, 0.765, 0.673, and 0.782 in GSE35640, VanAllen15, PUCH, Liu19, and Gide19 cohorts, respectively (Figure 7E7I), and was the highest in these cohorts except for GSE35640, suggesting a stable and robust performance in various independent datasets.

Comparison between ICDscore and other forms of cell death-based signatures. (A–D) The performance of ICDscore was compared with other forms of cell death-based signatures in predicting prognosis in the GSE65904 (A), TCGA-SKCM (B), GSE54467 (C), and GSE22153 (D) datasets. Statistic tests: two-sided z-score test. Data were presented as mean ± 95% confidence interval [CI]. (E–I) The performance of ICDscore was compared with other forms of cell death-based signatures in predicting immunotherapy efficiency in GSE35640 (E), VanAllen15 (F), PUCH (G), Liu19 (H), and Gide19 (I) cohorts. Ns, not significant; *p

Figure 7. Comparison between ICDscore and other forms of cell death-based signatures. (AD) The performance of ICDscore was compared with other forms of cell death-based signatures in predicting prognosis in the GSE65904 (A), TCGA-SKCM (B), GSE54467 (C), and GSE22153 (D) datasets. Statistic tests: two-sided z-score test. Data were presented as mean ± 95% confidence interval [CI]. (EI) The performance of ICDscore was compared with other forms of cell death-based signatures in predicting immunotherapy efficiency in GSE35640 (E), VanAllen15 (F), PUCH (G), Liu19 (H), and Gide19 (I) cohorts. Ns, not significant; *p < 0.05; **, p < 0.01; ***, p < 0.001; ****, p < 0.0001.

In addition to cell death-based signatures, many other gene expression signatures, for example, hypoxia- or CD8 T cell-related signatures, had been constructed to predict the prognosis of melanoma patients [24, 34]. By searching the PubMed website with the following query method “((signature[Title]) OR (classifier[Title])) AND (melanoma[Title])”, 52 mRNA signatures were ultimately enrolled for further comparison (Supplementary Table 4). These signatures were associated with DNA methylation, RNA methylation, oxidative stress, immune cell infiltration, metastasis, ferroptosis, proptosis or other biological process. Univariate Cox analyses revealed that only ICDscore, tumor immune-relevant (TIR) signature (Mei_2021) [35] and immune-associated genes (IAGs) signature (Song_2021) [36] had significant prognostic relevance across the four datasets (Figure 8A). The C-index of ICDscore displayed the best performance in GSE65904 and GSE54467 (Figure 8B, 8D), and ranked the fifth in the GSE22153 cohort (Figure 8E). In the TCGA-SKCM cohort, the C-index of ICDscore was not the highest, but was greater than 0.6 (Figure 8C). Besides, ICDscore had the highest average C-index (0.655) in all the four datasets (Figure 8F).

Comparison between ICDscore and other published signatures. (A) Univariate Cox regression analysis of ICDscore and 52 published signatures in GSE65904, TCGA-SKCM, GSE54467, and GSE22153 datasets. (B–E) C-index analyses of ICDscore and 52 published signatures in GSE65904 (B), TCGA-SKCM (C), GSE54467 (D), and GSE22153 (E) datasets. Statistic tests: two-sided z-score test. Data are presented as mean ± 95% confidence interval [CI]. (F) The average C-index of ICDscore and 52 published signatures across all the 4 datasets. Ns, not significant; *p

Figure 8. Comparison between ICDscore and other published signatures. (A) Univariate Cox regression analysis of ICDscore and 52 published signatures in GSE65904, TCGA-SKCM, GSE54467, and GSE22153 datasets. (BE) C-index analyses of ICDscore and 52 published signatures in GSE65904 (B), TCGA-SKCM (C), GSE54467 (D), and GSE22153 (E) datasets. Statistic tests: two-sided z-score test. Data are presented as mean ± 95% confidence interval [CI]. (F) The average C-index of ICDscore and 52 published signatures across all the 4 datasets. Ns, not significant; *p < 0.05; **, p < 0.01; ***, p < 0.001; ****, p < 0.0001.

Discussion

The capability of ICD to drive anti-cancer immune response provokes researchers’ interest in developing ICD-related risk model that predicts prognosis and response to immunotherapy in tumors [37, 38]. In both these two studies, ICD-associated subtyping systems were helpful for the development of precision immunotherapy. According to several published clinical studies, the ORR of previously untreated patients with metastatic melanoma is from about 10% to 57.6%, suggesting most of melanoma patients are insensitive to the ICIs [16, 18]. Given the tight association between ICD and immune response, ICD-related risk score might be helpful to identify melanoma patients potentially benefiting from immunotherapy.

In this work, two unsupervised clustering methods both showed that melanoma patients could be classified into two subgroups based on the expression of 34 ICD-related genes (Figure 2B2E). To minimize the potential bias of methodology, the clustering results of two methods were intersected, and melanoma patients grouped into C1 cluster by both methods were renamed as ‘Cluster A’, whereas those grouped into C2 cluster by both methods were renamed as ‘Cluster B’ (Figure 2H). The patients in the Cluster A had significantly elevated expression of most innate or adaptive immune system-related genes (Figure 2I). Consistently, patients in the Cluster A had a significant enrichment in most immune-related pathways like activation of immune response (Figure 3C) and was marked by dramatically higher infiltration of CD8 T cells (Figure 3G). Once ICD was triggered by anti-cancer therapy, the release of many damage-associated molecular patterns (DAMPs), such as cell surface-exposed calreticulin (CALR), extracellular adenosine triphosphate (ATP), could drive the activation and maturation of innate and subsequent adaptive immune cells [39]. Thus, it would be reasonable to propose that a higher expression of innate or adaptive immune system-related genes and infiltration of CD8 T cells might facilitate the occurrence of ICD and cause the amplification of anti-tumor effects induced by immune cells like cytotoxic T cells [37, 40, 41].

Further, the DEGs between Cluster A and B were input into a LASSO regression model, and 4 crucial genes were generated for the construction of ICDscore (Figure 4A4C). All the 4 genes participate in immune activities. GBP2 is a member of the p65 guanine-binding protein (GBP) family, which includes interferon-induced large GTPase that exhibits antiviral activity through the innate immune response [42, 43]. Lysozymes, encoded by LYZ, has antibacterial activity and constitutes the first line of defense [44, 45]. CST7 is predominantly expressed in immune cells and tends to attenuate the granule-dependent cytotoxicity of natural killer (NK) cells and T cells [4648]. SIRPG is expressed by T cells and interacts with CD47 on the surface of the cell, and modulates immune responses [49, 50]. Prognostic analyses in the testing dataset and three independent validating datasets revealed that high ICDscore correlated with significantly shorter OS and PFS (Figure 4D4G and Supplementary Figure 2B, 2C). In addition, ICDscore could serve as an independent prognostic factor (Figure 4I, 4J and Supplementary Table 5). Notably, ICDscore also had superior performance than other cell death related signatures in predicting prognosis of melanoma patients (Figure 7A7D). Moreover, we retrieved another 52 published mRNA signatures correlating with various biological activities. Univariate Cox regression displayed that only ICDscore and two other signatures maintained prognostic significance across all the four cohorts, suggesting most of the signatures have not been thoroughly validated (Figure 8A). It should be pointed out that many models had good performance in the training dataset but relatively poor performance in external validating datasets (Xu_2021-2 or Song_2020, for example) [51, 52]. Compared with other signatures, ICDscore exhibited stable performance across multiple cohorts and its average C-index was the highest, suggesting a better capability in predicting prognosis of melanoma patients (Figure 8B).

The value of ICDscore in guiding anti-cancer therapy was also investigated. On one hand, melanoma patients with low ICDscore were marked with significantly high infiltration of immune cells, particularly CD8 T cells (Figure 5B, 5D), suggesting an “immune-hot” phenotype (Figure 5E, 5F) [53]. These patients also exhibited high expression of ICI-related genes, such as PD-L1 (CD274) and CTLA4 (Figure 5G), whose high expression is generally associated with response to immunotherapy [54]. In addition, patients in the ICDscore-low subgroup also had higher TMB (Figure 5H), which could increase the production of mutation-derived neoantigens and contribute to activation of cytotoxic T cells [55]. In 6 independent cohorts in which melanoma patients receiving immunotherapy, patients who responded to immunotherapy displayed lower ICDscore (Figure 6A). Similarly, patients in the low-ICDscore subgroup generally had higher ORR to immunotherapy (Figure 6B) and longer OS (Figure 6C6F). The AUC of ICDscore in predicting response to immunotherapy was from 0.673 to 0.782, and was higher than that of other cell death related signatures (Figure 7E7I), suggesting a better performance than these signatures. On the other hand, ICDscore exhibited strong positive correlation with multiple targeted drugs like rapamycin and sunitinib (Figure 6G6N and Supplementary Table 3). Although some preclinical models or clinical trials had been conducted to evaluate the benefit of some of these drugs in melanoma [5658], the results are far from satisfactory in clinical practice. The ICDscore might be helpful in selecting suitable patients when conducting clinical trials.

At last, the limitation of this work should be pointed out that the ICDscore was developed and validated in a retrospective method; thus, evidence from a well-designed perspective analysis is required before the model be applied in clinics.

Conclusions

In conclusion, we constructed a stable and powerful ICD-related signature for evaluating the prognosis and benefits of immunotherapy, based on the integration of a set of bioinformatics tools. The ICDscore showed certain superiority than other mRNA signatures and served as a promising tool to guide decision-making and surveillance for individual melanoma patients.

Author Contributions

The study was designed by Guoyin Li, Mingsheng Yang and Changjing Wu. Data analysis was carried out by Guoyin Li, Huina Zhang, Jin Zhao, Qiongwen Liu and Jinke Jiao. Bioinformatics analysis was conducted by Guoyin Li. The manuscript was drafted by Guoyin Li, Mingsheng Yang and Changjing Wu, and was revised by all authors before the final version was approved to be published.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Funding

This work was supported by the National Natural Science Foundation of China (81903031), China Postdoctoral Science Foundation (2020M682334), Henan Postdoctoral Foundation (202003002), Fund of Key Laboratory of Modern Teaching Technology, Ministry of Education, China (NO.SYSK202107).

References

  • 1. Salaün H, de Koning L, Saint-Ghislain M, Servois V, Ramtohul T, Garcia A, Matet A, Cassoux N, Mariani P, Piperno-Neumann S, Rodrigues M. Nivolumab plus ipilimumab in metastatic uveal melanoma: a real-life, retrospective cohort of 47 patients. Oncoimmunology. 2022; 11:2116845. https://doi.org/10.1080/2162402X.2022.2116845 [PubMed]
  • 2. Schachter J, Ribas A, Long GV, Arance A, Grob JJ, Mortier L, Daud A, Carlino MS, McNeil C, Lotem M, Larkin J, Lorigan P, Neyns B, et al. Pembrolizumab versus ipilimumab for advanced melanoma: final overall survival results of a multicentre, randomised, open-label phase 3 study (KEYNOTE-006). Lancet. 2017; 390:1853–62. https://doi.org/10.1016/S0140-6736(17)31601-X [PubMed]
  • 3. Ahluwalia P, Kolhe R, Gahlay GK. The clinical relevance of gene expression based prognostic signatures in colorectal cancer. Biochim Biophys Acta Rev Cancer. 2021; 1875:188513. https://doi.org/10.1016/j.bbcan.2021.188513 [PubMed]
  • 4. Sparano JA, Gray RJ, Makower DF, Pritchard KI, Albain KS, Hayes DF, Geyer CE Jr, Dees EC, Perez EA, Olson JA Jr, Zujewski J, Lively T, Badve SS, et al. Prospective Validation of a 21-Gene Expression Assay in Breast Cancer. N Engl J Med. 2015; 373:2005–14. https://doi.org/10.1056/NEJMoa1510764 [PubMed]
  • 5. Paik S, Shak S, Tang G, Kim C, Baker J, Cronin M, Baehner FL, Walker MG, Watson D, Park T, Hiller W, Fisher ER, Wickerham DL, et al. A multigene assay to predict recurrence of tamoxifen-treated, node-negative breast cancer. N Engl J Med. 2004; 351:2817–26. https://doi.org/10.1056/NEJMoa041588 [PubMed]
  • 6. Peng F, Liao M, Qin R, Zhu S, Peng C, Fu L, Chen Y, Han B. Regulated cell death (RCD) in cancer: key pathways and targeted therapies. Signal Transduct Target Ther. 2022; 7:286. https://doi.org/10.1038/s41392-022-01110-y [PubMed]
  • 7. Kroemer G, Galassi C, Zitvogel L, Galluzzi L. Immunogenic cell stress and death. Nat Immunol. 2022; 23:487–500. https://doi.org/10.1038/s41590-022-01132-2 [PubMed]
  • 8. Li Z, Chu Z, Yang J, Qian H, Xu J, Chen B, Tian T, Chen H, Xu Y, Wang F. Immunogenic Cell Death Augmented by Manganese Zinc Sulfide Nanoparticles for Metastatic Melanoma Immunotherapy. ACS Nano. 2022; 16:15471–83. https://doi.org/10.1021/acsnano.2c08013 [PubMed]
  • 9. Zhang Q, Li S, Ren J, He X, Shi H, Zhang F, Li H, Tong R. ROS-triggered nanoinducer based on dermatan sulfate enhances immunogenic cell death in melanoma. J Control Release. 2022; 348:22–33. https://doi.org/10.1016/j.jconrel.2022.04.026 [PubMed]
  • 10. Ping S, Wang S, Zhao Y, He J, Li G, Li D, Wei Z, Chen J. Identification and validation of a ferroptosis-related gene signature for predicting survival in skin cutaneous melanoma. Cancer Med. 2022; 11:3529–41. https://doi.org/10.1002/cam4.4706 [PubMed]
  • 11. Fei H, Chen X. Establishment and validation of an autophagy-related prognostic signature for survival predicting in cutaneous melanoma. Am J Cancer Res. 2021; 11:5979–91. [PubMed]
  • 12. Ju A, Tang J, Chen S, Fu Y, Luo Y. Pyroptosis-Related Gene Signatures Can Robustly Diagnose Skin Cutaneous Melanoma and Predict the Prognosis. Front Oncol. 2021; 11:709077. https://doi.org/10.3389/fonc.2021.709077 [PubMed]
  • 13. Niu Z, Wang X, Xu Y, Li Y, Gong X, Zeng Q, Zhang B, Xi J, Pei X, Yue W, Han Y. Development and Validation of a Novel Survival Model for Cutaneous Melanoma Based on Necroptosis-Related Genes. Front Oncol. 2022; 12:852803. https://doi.org/10.3389/fonc.2022.852803 [PubMed]
  • 14. Tucci M, Passarelli A, Mannavola F, Felici C, Stucci LS, Cives M, Silvestris F. Immune System Evasion as Hallmark of Melanoma Progression: The Role of Dendritic Cells. Front Oncol. 2019; 9:1148. https://doi.org/10.3389/fonc.2019.01148 [PubMed]
  • 15. Han IH, Jeong C, Yang J, Park SH, Hwang DS, Bae H. Therapeutic Effect of Melittin-dKLA Targeting Tumor-Associated Macrophages in Melanoma. Int J Mol Sci. 2022; 23:3094. https://doi.org/10.3390/ijms23063094 [PubMed]
  • 16. Robert C, Schachter J, Long GV, Arance A, Grob JJ, Mortier L, Daud A, Carlino MS, McNeil C, Lotem M, Larkin J, Lorigan P, Neyns B, et al, and KEYNOTE-006 investigators. Pembrolizumab versus Ipilimumab in Advanced Melanoma. N Engl J Med. 2015; 372:2521–32. https://doi.org/10.1056/NEJMoa1503093 [PubMed]
  • 17. Hodi FS, O’Day SJ, McDermott DF, Weber RW, Sosman JA, Haanen JB, Gonzalez R, Robert C, Schadendorf D, Hassel JC, Akerley W, van den Eertwegh AJ, Lutzky J, et al. Improved survival with ipilimumab in patients with metastatic melanoma. N Engl J Med. 2010; 363:711–23. https://doi.org/10.1056/NEJMoa1003466 [PubMed]
  • 18. Larkin J, Chiarion-Sileni V, Gonzalez R, Grob JJ, Cowey CL, Lao CD, Schadendorf D, Dummer R, Smylie M, Rutkowski P, Ferrucci PF, Hill A, Wagstaff J, et al. Combined Nivolumab and Ipilimumab or Monotherapy in Untreated Melanoma. N Engl J Med. 2015; 373:23–34. https://doi.org/10.1056/NEJMoa1504030 [PubMed]
  • 19. Garg AD, De Ruysscher D, Agostinis P. Immunological metagene signatures derived from immunogenic cancer cell death associate with improved survival of patients with lung, breast or ovarian malignancies: A large-scale meta-analysis. Oncoimmunology. 2015; 5:e1069938. https://doi.org/10.1080/2162402X.2015.1069938 [PubMed]
  • 20. Cui C, Xu C, Yang W, Chi Z, Sheng X, Si L, Xie Y, Yu J, Wang S, Yu R, Guo J, Kong Y. Ratio of the interferon-γ signature to the immunosuppression signature predicts anti-PD-1 therapy response in melanoma. NPJ Genom Med. 2021; 6:7. https://doi.org/10.1038/s41525-021-00169-w [PubMed]
  • 21. Zhu Z, Li G, Li Z, Wu Y, Yang Y, Wang M, Zhang H, Qu H, Song Z, He Y. Core immune cell infiltration signatures identify molecular subtypes and promote precise checkpoint immunotherapy in cutaneous melanoma. Front Immunol. 2022; 13:914612. https://doi.org/10.3389/fimmu.2022.914612 [PubMed]
  • 22. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. 2010; 26:1572–3. https://doi.org/10.1093/bioinformatics/btq170 [PubMed]
  • 23. Gaujoux R, Seoighe C. A flexible R package for nonnegative matrix factorization. BMC Bioinformatics. 2010; 11:367. https://doi.org/10.1186/1471-2105-11-367 [PubMed]
  • 24. Yuan Y, Zhu Z, Lan Y, Duan S, Zhu Z, Zhang X, Li G, Qu H, Feng Y, Cai H, Song Z. Development and Validation of a CD8+ T Cell Infiltration-Related Signature for Melanoma Patients. Front Immunol. 2021; 12:659444. https://doi.org/10.3389/fimmu.2021.659444 [PubMed]
  • 25. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, Treviño V, Shen H, Laird PW, Levine DA, Carter SL, Getz G, Stemke-Hale K, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013; 4:2612. https://doi.org/10.1038/ncomms3612 [PubMed]
  • 26. Jun Z. GseaVis. An Implement R Package to Visualize GSEA Results. 2022.
  • 27. Kist M, Vucic D. Cell death pathways: intricate connections and disease implications. EMBO J. 2021; 40:e106700. https://doi.org/10.15252/embj.2020106700 [PubMed]
  • 28. Network T. The Immune Landscape of Cancer. 2018.
  • 29. Cancer Genome Atlas Network. Genomic Classification of Cutaneous Melanoma. Cell. 2015; 161:1681–96. https://doi.org/10.1016/j.cell.2015.05.044 [PubMed]
  • 30. Bagaev A, Kotlov N, Nomie K, Svekolkin V, Gafurov A, Isaeva O, Osokin N, Kozlov I, Frenkel F, Gancharova O, Almog N, Tsiper M, Ataullakhanov R, Fowler N. Conserved pan-cancer microenvironment subtypes predict response to immunotherapy. Cancer Cell. 2021; 39:845–65.e7. https://doi.org/10.1016/j.ccell.2021.04.014 [PubMed]
  • 31. Chen D, Zhang Y, Wang W, Chen H, Ling T, Yang R, Wang Y, Duan C, Liu Y, Guo X, Fang L, Liu W, Liu X, et al. Identification and Characterization of Robust Hepatocellular Carcinoma Prognostic Subtypes Based on an Integrative Metabolite-Protein Interaction Network. Adv Sci (Weinh). 2021; 8:e2100311. https://doi.org/10.1002/advs.202100311 [PubMed]
  • 32. Zheng C, Zheng L, Yoo JK, Guo H, Zhang Y, Guo X, Kang B, Hu R, Huang JY, Zhang Q, Liu Z, Dong M, Hu X, et al. Landscape of Infiltrating T Cells in Liver Cancer Revealed by Single-Cell Sequencing. Cell. 2017; 169:1342–56.e16. https://doi.org/10.1016/j.cell.2017.05.035 [PubMed]
  • 33. Hodi FS, Wolchok JD, Schadendorf D, Larkin J, Long GV, Qian X, Saci A, Young TC, Srinivasan S, Chang H, Tang H, Wind-Rotolo M, Rizzo JI, et al. TMB and Inflammatory Gene Expression Associated with Clinical Outcomes following Immunotherapy in Advanced Melanoma. Cancer Immunol Res. 2021; 9:1202–13. https://doi.org/10.1158/2326-6066.CIR-20-0983 [PubMed]
  • 34. Shou Y, Yang L, Yang Y, Zhu X, Li F, Xu J. Determination of hypoxia signature to predict prognosis and the tumor immune microenvironment in melanoma. Mol Omics. 2021; 17:307–16. https://doi.org/10.1039/d0mo00159g [PubMed]
  • 35. Mei Y, Chen MM, Liang H, Ma L. A four-gene signature predicts survival and anti-CTLA4 immunotherapeutic responses based on immune classification of melanoma. Commun Biol. 2021; 4:383. https://doi.org/10.1038/s42003-021-01911-x [PubMed]
  • 36. Song LB, Luan JC, Zhang QJ, Chen L, Wang HY, Cao XC, Song NH, Lu Y. The Identification and Validation of a Robust Immune-Associated Gene Signature in Cutaneous Melanoma. J Immunol Res. 2021; 2021:6686284. https://doi.org/10.1155/2021/6686284 [PubMed]
  • 37. Wang X, Wu S, Liu F, Ke D, Wang X, Pan D, Xu W, Zhou L, He W. An Immunogenic Cell Death-Related Classification Predicts Prognosis and Response to Immunotherapy in Head and Neck Squamous Cell Carcinoma. Front Immunol. 2021; 12:781466. https://doi.org/10.3389/fimmu.2021.781466 [PubMed]
  • 38. Xu M, Lu JH, Zhong YZ, Jiang J, Shen YZ, Su JY, Lin SY. Immunogenic Cell Death-Relevant Damage-Associated Molecular Patterns and Sensing Receptors in Triple-Negative Breast Cancer Molecular Subtypes and Implications for Immunotherapy. Front Oncol. 2022; 12:870914. https://doi.org/10.3389/fonc.2022.870914 [PubMed]
  • 39. Serrano-Del Valle A, Anel A, Naval J, Marzo I. Immunogenic Cell Death and Immunotherapy of Multiple Myeloma. Front Cell Dev Biol. 2019; 7:50. https://doi.org/10.3389/fcell.2019.00050 [PubMed]
  • 40. Qiu X, Qu Y, Guo B, Zheng H, Meng F, Zhong Z. Micellar paclitaxel boosts ICD and chemo-immunotherapy of metastatic triple negative breast cancer. J Control Release. 2022; 341:498–510. https://doi.org/10.1016/j.jconrel.2021.12.002 [PubMed]
  • 41. Krombach J, Hennel R, Brix N, Orth M, Schoetz U, Ernst A, Schuster J, Zuchtriegel G, Reichel CA, Bierschenk S, Sperandio M, Vogl T, Unkel S, et al. Priming anti-tumor immunity by radiotherapy: Dying tumor cell-derived DAMPs trigger endothelial cell activation and recruitment of myeloid cells. Oncoimmunology. 2018; 8:e1523097. https://doi.org/10.1080/2162402X.2018.1523097 [PubMed]
  • 42. MacMicking JD. IFN-inducible GTPases and immunity to intracellular pathogens. Trends Immunol. 2004; 25:601–9. https://doi.org/10.1016/j.it.2004.08.010 [PubMed]
  • 43. Zhang S, Chen K, Zhao Z, Zhang X, Xu L, Liu T, Yu S. Lower Expression of GBP2 Associated With Less Immune Cell Infiltration and Poor Prognosis in Skin Cutaneous Melanoma (SKCM). J Immunother. 2022. https://doi.org/10.1097/CJI.0000000000000421 [PubMed]
  • 44. Yadav S, Surolia A. Lysozyme elicits pain during nerve injury by neuronal Toll-like receptor 4 activation and has therapeutic potential in neuropathic pain. Sci Transl Med. 2019; 11:eaav4176. https://doi.org/10.1126/scitranslmed.aav4176 [PubMed]
  • 45. Lelouard H, Henri S, De Bovis B, Mugnier B, Chollat-Namy A, Malissen B, Méresse S, Gorvel JP. Pathogenic bacteria and dead cells are internalized by a unique subset of Peyer’s patch dendritic cells that express lysozyme. Gastroenterology. 2010; 138:173–84.e1. https://doi.org/10.1053/j.gastro.2009.09.051 [PubMed]
  • 46. Kos J, Nanut MP, Prunk M, Sabotič J, Dautović E, Jewett A. Cystatin F as a regulator of immune cell cytotoxicity. Cancer Immunol Immunother. 2018; 67:1931–8. https://doi.org/10.1007/s00262-018-2165-5 [PubMed]
  • 47. Prunk M, Nanut MP, Sabotic J, Svajger U, Kos J. Increased cystatin F levels correlate with decreased cytotoxicity of cytotoxic T cells. Radiol Oncol. 2019; 53:57–68. https://doi.org/10.2478/raon-2019-0007 [PubMed]
  • 48. Maher K, Konjar S, Watts C, Turk B, Kopitar-Jerala N. Cystatin F regulates proteinase activity in IL-2-activated natural killer cells. Protein Pept Lett. 2014; 21:957–65. https://doi.org/10.2174/0929866521666140403124146 [PubMed]
  • 49. Nettleship JE, Ren J, Scott DJ, Rahman N, Hatherley D, Zhao Y, Stuart DI, Barclay AN, Owens RJ. Crystal structure of signal regulatory protein gamma (SIRPγ) in complex with an antibody Fab fragment. BMC Struct Biol. 2013; 13:13. https://doi.org/10.1186/1472-6807-13-13 [PubMed]
  • 50. Sinha S, Borcherding N, Renavikar PS, Crawford MP, Tsalikian E, Tansey M, Shivapour ET, Bittner F, Kamholz J, Olalde H, Gibson E, Karandikar NJ. An autoimmune disease risk SNP, rs2281808, in SIRPG is associated with reduced expression of SIRPγ and heightened effector state in human CD8 T-cells. Sci Rep. 2018; 8:15440. https://doi.org/10.1038/s41598-018-33901-1 [PubMed]
  • 51. Xu Z, Xie Y, Mao Y, Huang J, Mei X, Song J, Sun Y, Yao Z, Shi W. Ferroptosis-Related Gene Signature Predicts the Prognosis of Skin Cutaneous Melanoma and Response to Immunotherapy. Front Genet. 2021; 12:758981. https://doi.org/10.3389/fgene.2021.758981 [PubMed]
  • 52. Song LB, Zhang QJ, Hou XY, Xiu YY, Chen L, Song NH, Lu Y. A twelve-gene signature for survival prediction in malignant melanoma patients. Ann Transl Med. 2020; 8:312. https://doi.org/10.21037/atm.2020.02.132 [PubMed]
  • 53. Galon J, Bruni D. Approaches to treat immune hot, altered and cold tumours with combination immunotherapies. Nat Rev Drug Discov. 2019; 18:197–218. https://doi.org/10.1038/s41573-018-0007-y [PubMed]
  • 54. Song Y, Li Z, Xue W, Zhang M. Predictive biomarkers for PD-1 and PD-L1 immune checkpoint blockade therapy. Immunotherapy. 2019; 11:515–29. https://doi.org/10.2217/imt-2018-0173 [PubMed]
  • 55. McGranahan N, Furness AJ, Rosenthal R, Ramskov S, Lyngaa R, Saini SK, Jamal-Hanjani M, Wilson GA, Birkbak NJ, Hiley CT, Watkins TB, Shafi S, Murugaesu N, et al. Clonal neoantigens elicit T cell immunoreactivity and sensitivity to immune checkpoint blockade. Science. 2016; 351:1463–9. https://doi.org/10.1126/science.aaf1490 [PubMed]
  • 56. Decoster L, Vande Broek I, Neyns B, Majois F, Baurain JF, Rottey S, Rorive A, Anckaert E, De Mey J, De Brakeleer S, De Grève J. Biomarker Analysis in a Phase II Study of Sunitinib in Patients with Advanced Melanoma. Anticancer Res. 2015; 35:6893–9. [PubMed]
  • 57. Minor DR, Kashani-Sabet M, Garrido M, O’Day SJ, Hamid O, Bastian BC. Sunitinib therapy for melanoma patients with KIT mutations. Clin Cancer Res. 2012; 18:1457–63. https://doi.org/10.1158/1078-0432.CCR-11-1987 [PubMed]
  • 58. Kenessey I, Kramer Z, István L, Cserepes MT, Garay T, Hegedűs B, Dobos J, Tímár J, Tóvári J. Inhibition of epidermal growth factor receptor improves antitumor efficacy of vemurafenib in BRAF-mutant human melanoma in preclinical model. Melanoma Res. 2018; 28:536–46. https://doi.org/10.1097/CMR.0000000000000488 [PubMed]