Introduction

Melanoma, a highly malignant form of skin cancer, originates from abnormal melanocytes [1]. According to the Global Cancer Statistics 2022, approximately 330,000 new cases of melanoma were reported worldwide, with around 59,000 deaths resulting from the disease. [2]. Patients diagnosed with early-stage melanoma often have a high probability of cure through surgical intervention, with a five-year survival rate exceeding 90%. However, late-stage or metastatic melanoma typically presents a poor prognosis and the therapeutic options are considerably limited [2–4]. Currently, immune checkpoint blockade (ICB) therapies have achieved transformative breakthroughs in the treatment of various cancers, with notable success in treating metastatic melanoma. Over recent decades, the development of immune checkpoint inhibitors that target cytotoxic T lymphocyte-associated antigen 4 (CTLA-4), programmed death receptor 1 (PD-1), and programmed death-ligand 1 (PD-L1) has delivered promising outcomes for patients suffering from malignant melanoma [5]. Despite improvements in the prognosis of patients with metastatic melanoma through ICB therapy, a significant subset remains unresponsive to this treatment, deriving no benefit from it. Approximately half of these patients experience relapse and metastasis due to resistance to the therapy [6]. Consequently, it is crucial to delve deeper into the factors that influence the efficacy of ICB therapy. There is a pressing need for more accurate predictive biomarkers to determine the effectiveness of ICB and to devise more precise and effective treatment strategies for patients.

Currently, while PD-L1 expression, tumor mutational burden (TMB), and microsatellite instability-high (MSI-H) hold promise as predictive biomarkers for immunotherapy response in melanoma, they also have limitations related to consistency, heterogeneity, cut-off selection, prevalence, and response variability [7–9]. Additionally, research has identified novel biomarkers such as the interferon-γ signaling pathway genes and IGF1R mutation that can also be used to predict the effectiveness of immunotherapy in melanoma [10–11]. However, their predictive accuracy and clinical utility need further validation.

Immunogenic cell death (ICD) represents a type of tumor cell death induced by stressors like specific chemotherapy agents, oncolytic viruses, photodynamic therapy, and radiation [12–13]. The initiation of ICD leads to enhanced immunogenicity in tumor cells and the release of damage-associated molecular patterns (DAMPs), which in turn activate T cell-mediated specific anti-tumor immune responses [12–13]. The discovery of ICD has paved new pathways for cancer therapy, emphasizing that inducing ICD can amplify the immunogenic properties of tumor treatments. Notably, two potent ICD inducers, lurbinectedin and belantamab mafodotin, have received FDA approval for use in anti-tumor therapies in humans [14, 15]. Moreover, studies have identified DAMPs as having predictive value in assessing responses to immune therapy among cancer patients [16]. Consequently, investigating ICD-related biomarkers in melanoma and developing a predictive model for immune therapy efficacy could significantly enhance the accuracy and effectiveness of immune therapies for melanoma, presenting substantial clinical relevance. However, whether ICD signatures can predict the immune therapy response in on-treatment metastatic melanoma specimens is currently unclear. In this study, we utilized the elastic net algorithm to identify 9 hub ICD signatures as effective indicators for predicting the efficacy of ICB therapy. These markers include Alkaliptosis, Apoptosis, Cuproptosis, Entosis, ICD, Lysosome-dependent cell death (LDCD), Necroptosis, Oxeiptosis, and Parthanatos. Each form of cell death has distinct mechanisms but often converges on immune modulation. For instance, Alkaliptosis, Cuproptosis, Oxeiptosis, and Parthanatos are stress-induced cell death pathways, each involving distinct mechanisms such as oxidative damage or metal dysregulation, which can subsequently lead to immune activation. [17–21]. Apoptosis and Necroptosis trigger immune responses via caspase activation and inflammatory signaling, respectively [22–24]. Entosis and LDCD modulate immune responses primarily through promoting cell clearance and releasing lysosomal enzymes, which may alter tumor immunogenicity [25, 26]. Understanding the mechanisms of these forms of cell death is crucial for developing effective immunotherapies against melanoma.

Methods

Data collection and patients selection

To assess the effectiveness of ICB therapy in patients with metastatic melanoma, we collated data from several sources: Riaz et al. cohort (GEO: GSE120575), Gide et al. cohort (BioProject: PRJEB23709), MGH cohort (GEO: GSE115821 and GSE168204), and Abril-Rodriguez et al. cohort (dbGaP: phs001919.v1.p1). We included the on-treatment samples and excluded those sequencing data were incomplete or lacked sufficient treatment response assessments. All RNA gene raw counts were normalized to transcripts per kilobase million (TPM) expression. A robust multi-array quantile method was then employed to normalize the data from the different arrays. Utilizing the RECIST criteria [27], we classified the patients into two groups: responders (R) and nonresponders (NR). R were defined as those achieving complete response (CR), partial response (PR), or maintaining stable disease (SD) with a progression-free survival (PFS) exceeding 180 days. NR were identified as those with a PFS under 180 days and exhibiting disease progression (PD).

Calculation of signature scores

In this study, we leverage the 18 ICD gene signatures (Supplementary Table S1) for calculating the ICD scores for each specimen utilizing the single sample gene set enrichment analysis (ssGSEA) approach. Subsequently, the ‘glmnet’ R package was deployed to apply an Elastic-Net penalized logistic regression modelel for identifying ICD signatures associated with the efficacy of ICB therapy [28]. To mitigate overfitting, a three-fold cross-validation was implemented on the training set. We addressed the challenge of imbalanced classification through a cost-sensitive algorithm. The effectiveness of the predictive model was evaluated by generating receiver operating characteristic (ROC) curves. The predictive accuracy of the ICD signature scores in forecasting ICB therapy outcomes was quantified by the area under the curve (AUC). Using data from Riaz et al. cohort, the optimal cutoff was determined via Youden’s index method from the ROC analysis [29, 30]. Finally, we calculated the sample’s odd ratio according to the ICD signature score [31].

Statistical analysis

We employed the Wilcoxon one-sided rank-sum test to evaluate the statistical differences between R and NR within the ICB cohort. The samples were stratified into low- and high-signature score group according to the predetermined cutoff value. We assessed the survival disparities between these groups using Kaplan-Meier (KM) estimates and log-rank tests, and samples were separated into high and low group use median value of samples’ odd ratio as cutoff. Additionally, hazard ratios (HR) for all models were determined via Cox regression analysis. The verification of the proportional hazards assumption for the Cox model across all cohorts is presented in the Supplementary Fig. 1. All statistical analyses and graphics were performed using R software (version 4.2.1) and GraphPad Prism (version 8.0).

Result

Data collection and processing

We gathered RNA transcription data and clinical details from four datasets evaluating the responses to ICB therapy in patients with metastatic melanoma, including the Riaz et al. cohort, Gide et al. cohort, MGH cohort, and Abril-Rodriguez et al. cohort. As depicted in Fig. 1, the Riaz et al. cohort was chosen as the training set, incorporating 54 on-treatment samples for constructing an Elastic Net Regression (ELNR) model. The cohorts from Gide et al., MGH, and Abril-Rodriguez et al. were utilized as validation sets.

Fig. 1
figure 1

The flow chart of ICD signature score construction.

Full size image

We initially computed the ICD scores for each sample in the Riaz et al. training dataset using the ssGSEA algorithm, based on 18 distinct ICD gene signatures. The ENLR model was then employed to pinpoint the ICD signatures that offered the highest predictive accuracy. To prevent overfitting, three-fold cross-validation was applied, and a cost-sensitive algorithm was utilized to handle the issue of imbalanced classification (Fig. 2A and B). Ultimately, 9 hub ICD signatures were identified as the most effective predictors of the response to ICB therapy. These include Alkaliptosis, Apoptosis, Cuproptosis, Entosis, ICD, LDCD, Necroptosis, Oxeiptosis, and Paraptosis (Fig. 2C). Additionally, we leveraged the effective sizes as weighting factors and computed the weighted averages of the ssGSEA scores for the 9 ICD signatures, designating this metric as the ICD-based signature scores (ICDS) for on-treatment samples. Next, based on the average odds ratio derived from the ICDS, the patients in the training cohort were divided into high and low score groups.

Fig. 2
figure 2

(A, B). The model’s training parameter selection process to generate Riaz et al. on-treatment samples for the ICD-based signatures. To avoid overfitting, 3-fold cross-validation was performed with the parameter setting as ʺtype.measure = auc, family = ‘binomial’.ʺ (C). ICD-based signatures consists of 9 selected ICD associated with the effect sizes (variable weights) from elastic net penalized logistic regression model.

Full size image

ICD-based signatures for on-treatment samples

In the Riaz et al. training dataset, we generated a heatmap for each on-treatment NR and R to illustrate the ssGSEA scores across the 9 ICD signatures (Fig. 3A). In the training dataset, patients with high ICDS demonstrated a notably higher response rate to ICB therapy compared to the low ICDS group, with rates of 54.2% and 20.8%, respectively (Fig. 3B). A comparison between the R group and NR group showed a significantly higher ICDS in the R group (p < 0.001) (Fig. 3C). ROC analysis indicated that the AUC for differentiating between R and NR based on ICDS was 0.81 in the training set (Fig. 3D). Furthermore, survival analysis revealed that patients with higher ICDS exhibited longer progression-free survival (PFS) and overall survival (OS) than those with lower ICDS (Fig. 3E and F).

Fig. 3
figure 3

(A) Heatmap representing the single sample gene set enrichment analysis value of on-treatment NR and R in the Riaz et al. cohort. (B) The response rates of the low ICDS group and the high ICDS group for on-treatment samples from the Riaz et al. cohort. (C) Boxplot of ICDS for on-treatment samples from the Riaz et al. cohort. P values were computed via a one-sided Wilcoxon rank-sum test. (D) Receiver operating correlation curve and area under the curve of ICD signatures for on-treatment samples from the Riaz et al. cohort. (E, F) Kaplan–Meier curves of PFS and OS for on-treatment samples based on ICD signature scores for the Riaz et al. cohort.

Full size image

Additionally, we generated heatmaps for each on-treatment NR and R to display the ssGSEA scores for the 9 ICD signatures across three validation cohorts (Fig. 4A). Our findings revealed that, in these cohorts, R’s had significantly higher ICDS compared to NR’s (Fig. 4B). Furthermore, patients with high ICDS demonstrated a significantly higher response rate to ICB treatment than those with low ICDS (Fig. 4C). The AUCs for the Gide et al., AbrilRodriguez et al., and MGH et al. datasets were 0.83, 0.81, and 0.79, respectively (Fig. 4D), showcasing the model’s predictive capability. Finally, through survival analysis of the Gide et al. and MGH et al. cohorts, we observed that patients with high ICDS had longer OS and PFS, although no statistical difference in OS was noted between high and low ICDS patients in the Gide et al. cohort (Fig. 4E and F).

Fig. 4
figure 4

(A) Heatmap representing the single sample gene set enrichment analysis value of on-treatment NR and R in the MGH cohort, the Abril-Rodriguez et al. cohort and the Gide et al. cohort. (B) The response rates of the low ICDS group and the high ICDS group for on-treatment samples from the MGH cohort, the Abril-Rodriguez et al. cohort and the Gide et al. cohort. (C) Boxplot of ICDS for on-treatment samples from the MGH cohort, the Abril-Rodriguez et al. cohort and the Gide et al. cohort. P values were computed via a one-sided Wilcoxon rank-sum test. (D) Receiver operating correlation curve and area under the curve of ICD signatures for on-treatment samples from the MGH cohort, the Abril-Rodriguez et al. cohort and the Gide et al. cohort. (E, F) Kaplan–Meier curves of PFS and OS for on-treatment samples based on ICD signature scores for the MGH cohort and the Gide et al. cohort.

Full size image

ICDS was an independent prognostic factor

To determine whether the ICDS was an independent prognostic factor for predicting the prognosis of patients with metastatic melanoma, we conducted a multivariate Cox analysis in the Riaz et al. training dataset, including four variables: ICDS, tumor type, mutation type, and tumor stage. The results indicate that ICDS is an independent prognostic marker for predicting the PFS (HR = 0.15, p = 0.03) and OS (HR = 0.01, p = 0.01) of patients with metastatic melanoma (Fig. 5A and B). Similarly, in the Gide et al. validation set, we included variables such as ICD, age, and gender for analysis. In this validation set, ICDS independently predicts PFS (HR = 0.65, p = 0.03) (Fig. 5C). However, the HR for predicting OS using ICDS was 0.59, which did not reach statistical significance (p = 0.08) (Fig. 5D).

Fig. 5
figure 5

(A, B) Multi-variate regression analysis of PFS (A) and OS (B) in on-treatment samples patients from Riaz et al. cohort. (C, D) Multi-variate regression analysis of PFS (C) and OS (D) in on-treatment samples patients from Gide et al. cohort.

Full size image

Comparison of ICD signatures and published predictive signatures

We evaluated the predictive performance of the ICD signatures derived from on-treatment specimens. To this end, we further compared the ICD signatures scores against those of previously reported transcriptome-based predictive signatures, including the PDL1. Sig [32], C-X-C Motif Chemokine Ligand (CXCL). Sig [33], IRG. Sig [34], Interferon Gamma (IFNG). Sig [35], immuno-predictive score (IMPRES). Sig [36], Blood. Sig [37], innate anti-PD-1 resistance (IPRES). Sig [38], anti-CTLA4 resistance MAGE gene (CRAM). Sig [39], epithelial-to-mesenchymal transition (EMT). Sig [40], LRRC15 + carcinoma-associated fibroblasts (LRRC15. CAF). Sig [41]. ICD signatures achieved average AUCs of 0.81, which were superior to any other published signatures (AUCs ranging from 0.58 to 0.76) (Fig. 6).

Fig. 6
figure 6

Comparison for the predictive ICB treatment response performance of ICD signatures and published predictive signatures.

Full size image

Discussion

Immunotherapy represents a promising frontier in cancer treatment, offering new hope for patients [42–43]. Currently, many biomarkers are used to predict the response to immunotherapy in melanoma, such as PD-L1 expression, TMB, and MSI-H. However, they also have limitations such as consistency, heterogeneity and cut-off selection [7–9]. Therefore, the construction of predictive models for immunotherapy has significant clinical implications. Construction predictive models for immunotherapy not only enhances the accuracy of predicting individual patient responses, thereby facilitating more tailored treatment plans and improving treatment outcomes, but it also deepens our understanding of the biological characteristics of tumors and their interactions with the immune system [44–48]. This contributes significantly to the development of innovative therapeutic strategies for cancer and steers the direction of future scientific research. ICD plays a significant role in cancer immunotherapy. ICD can induce the production of immunostimulatory factors, such as ATP, IL-1β, type I interferons, CXCL10, and HMGB1 [49]. These molecules can regulate the composition and function of innate and adaptive immune cells within the tumor microenvironment. In vitro studies have shown that various types of antitumor treatments can upregulate stress-induced molecules in tumor cells, making them more susceptible to NK cell-mediated killing [49]. Recent studies have shown that stimulating ICD through nanotechnology and intelligent responsive systems can significantly enhance cancer immunotherapy [50]. For example, Academician Tang Benzhong’s team from the Chinese University of Hong Kong has developed a light-triggered self-accelerating nano-platform for multifunctional image-guided combination cancer immunotherapy, bringing great hope for cancer treatment [51]. Additionally, in recent years, many researchers have developed risk scoring models based on ICD signatures [52–54], which means that the expression levels of these genes can be used to predict tumor prognosis. However, there is currently no research on whether the ICD gene signature can predict the response and outcome of on-treatment metastatic melanoma patients to ICB therapy. In this study, we utilized 18 ICD gene signatures and employed an elastic net algorithm combined with four distinct datasets to establish a predictive model. This model aims to predict the response of metastatic melanoma patients to immunotherapy.

Through the ENLR algorithm, we ultimately identified 9 hub ICD signatures as the most effective predictors of the response to ICB therapy. These markers include Alkaliptosis, Apoptosis, Cuproptosis, Entosis, ICD, LDCD, Necroptosis, Oxeiptosis, and Parthanatos. Based on 9 ICD signatures, we analyzed on-treatment samples of metastatic melanoma across four independent cohorts. We successfully constructed and validated a clinical prediction model to predict the ICB response in patients with metastatic melanoma. Our data indicate that metastatic melanoma patients with higher ICDS are more likely to benefit from ICB therapy. Through ROC analysis, we found that both the training cohort and the validation cohort showed an AUC of around 0.8, demonstrating the effective predictive performance of this model. Additionally, our research indicates that patients with higher ICDS have better prognoses, and ICDS is an independent predictor for patients with metastatic melanoma.

In the field of predicting the efficacy of ICB therapy, an increasing number of biomarkers have been developed and validated. For instance, Du et al. developed pathway-based super signatures (PASS) to predict the response of metastatic melanoma patients to anti-PD1 ICB therapy [48]. Although our study and Du et al. study utilized similar datasets and analysis methods, we specifically highlights the role of ICD in predicting the efficacy of immune checkpoint blockade ICB therapy. Quite the contrary, Du et al. study did not study or mention about ICD signatures. The ICD gene signatures have not been previously described for predicting response of ICD therapy, particularly in melanoma samples undergoing treatment. By introducing the ICD gene signature, our research uniquely links the cell death pathway to the prediction of ICB efficacy, emphasizing the critical role of stress-induced cell death in anti-tumor immune responses. This novel focus on ICD enriches the understanding of ICB mechanisms and offers new avenues for enhancing tumor immunity by targeting cell death pathways. While the predictive performance of Du et al.‘s study is comparable to ours, our research exhibits great potential of ICD signatures for clinical translation. Moreover, we conducted a Cox proportional hazards analysis that Du et al. study didn’t consider, integrating the ICD score with clinical factors, and demonstrated that the ICD score is an independent prognostic marker. Our analysis also examined the response rates between high and low ICD score groups, an analysis not covered by Du et al.‘s study. Above analysis provides deeper insights into how ICD-related mechanisms correlate with patient outcomes, enhancing the potential for real-time treatment adjustments based on our predictive model.

When compared to previously published signatures [32–41], ICD signature showed more robust and stable predictive performance. It is now widely accepted that on-treatment tumor specimens may provide more valuable insights into the dynamic changes at the transcriptional level correlated with clinical response, resulting in a higher predictive score. To a certain extent, these advances might explain why ICD signature is superior to previously published signatures, which are derived from pretreatment tumor specimens. These findings have implications for treatment selection in current clinical practice because most studies have used pretreatment biopsies to construct prediction signature. Even though we have demonstrated the strength performance of ICD-ON signature, several limitations of this signature still exist, which remain to be further refined. First, owing to dataset access limitations, only four cohorts of metastatic melanoma were obtained and analyzed in the study. The current study included a relatively small sample size and was retrospective in nature. Further validation studies in independent cohorts are essential before it can be implemented for clinical practice. Secondly, scores induced by ssGSEA method are potentially relative and highly sensitive to different variables such as batch effect, sequencing method, library preparation, etc. Finally, further studies are warranted to test the predictive performance of this signature in other types of cancers and ICB therapies such as CAR-T or oncolytic viral therapies.

In summary, the study of ICD in the immunotherapy of metastatic melanoma has shown great potential. Not only does it play a crucial role in stimulating anti-tumor immune responses, but it can also predict ICB treatment responses and prognosis more accurately through the expression of ICD signatures, thereby providing a basis for personalized treatment.

Conclusion

By employing the ENLR algorithm, we identified 9 ICD signatures that are closely associated with the immune response to ICB therapy in metastatic melanoma samples undergoing treatment. These signatures consistently and accurately predicted on-treatment patients’ responses to ICB treatment. These ICD signatures not only forecasted the prognosis of ICB therapy but also provided insights for clinical decision-making in on-treatment patients with metastatic melanoma.