Genetic variation landscape of necroptosis genes in cutaneous melanoma

We explored the genetic variation of necroptosis-related genes in cutaneous melanoma from 3 aspects in the TCGA database: CNV, somatic mutations, and differential expression in the transcriptome. We found that a similar number of necroptosis-related genes gained copy number and CNV deletion (Fig. 1A). Genomic mutations were common in these genes with 376 (80.51%) of 467 patients having occurred genetic changes, in which BRAF (51%) had the most genetic alteration. CDKN2A and HDAC9 followed with 12% mutational frequency (Fig. 1B). The location of the 65 necroptosis-related genes in human chromosomes could be seen in Fig. 1C. The differential analysis in normal tissue and tumor tissue showed 14 necroptosis-related genes with significantly differential expression in cutaneous melanoma (Fig. 1D). The alteration and genetic variation of necroptosis-related genes act as an important part in regulating the happening, aggravation, and prognosis of cutaneous melanoma.

DEGs between normal and tumor tissue specimens

The 65 necroptosis-associated genes were comparatively assessed for their expression in the pooled GTEx and TCGA data involving 813 normal and 473 tumor tissue specimens, and 14 DEGs were identified (P < 0.05). Of these DEGs, 6 (EGFR, GATA3, DIABLO, CFLAR, RIPK3, and MAPK8) were repressed, while (ZBP1, PLK1, SPATA2, BCL2, ALK, FASLG, TNFRSF21, and LEF1) were upregulated in cancer cases (Fig. 1D).

Identification of necroptosis patterns based on DEGs

Based on the expression profile of 14 DEGs, NMF consensus clustering analysis was carried out. At a clustering variable (k) of 2, the magnitude of the cophenetic correlation coefficient begins to fall, and intragroup and intergroup correlations were high and low, respectively, suggesting the SKCM cohort clustered in two groups (C1 and C2; Supplementary Fig. 1). We marked them necroptosis cluster A (n = 171) and B (n = 283), respectively (Fig. 2A; Table S4). The Kaplan–Meier curves showed that OS in the necroptosis cluster A was significantly better than that in the necroptosis cluster B (Fig. 2B).

Identification of necroptosis patterns based on DEGs. (A) NMF consensus clustering for the k value was 2. (B) Kaplan–Meier curves of OS for patients with the two clusters. (C) The distinction in 22 kinds of immune infiltration cells between the two clusters.

Immune cells infiltration analysis in TME

To understand the characteristics of immune infiltration cells in different clusters, we conducted the distinction of 22 kinds of immune cell infiltration between the two clusters. Among the immune infiltration cells, a majority of immune cells were found infiltrated into the necroptosis cluster A, including activated memory CD4 + T cells, CD8 + T cells, plasma cells, activated NK cells, monocytes, M1 macrophages, helper follicular T cells, and regulatory T cells (Tregs) (Fig. 2C).

Prognostic gene-based model in TCGA cases and external validation of the risk signature

Univariate Cox proportional hazards regression analysis showed that 7 out of the 14 necroptosis-related genes had P < 0.05 in the TCGA cohort (Fig. 3A). To eliminate collinearity of the variables and avoid overfitting of the prognostic model, these 7 genes were undergone the LASSO regression analysis in the TCGA dataset (Fig. 3B and C). Subsequently, 6 candidate genes were identified for further multivariate Cox regression analysis. Finally, the signature was constructed according to the expression levels of 4 genes (ZBP1, PLK1, EGFR, TNFRSF21, and ACAT2) (Fig. 3D). Risk score = -− 0.421 × (ZBP1 expression value) + 0.242 × (PLK1 expression value) + 0.169 × (EGFR expression value) -0.095 × (TNFRSF21 expression value). Based on the obtained median score, patients were divided into the low- and high-risk subgroups. In the TCGA and validation cohorts, high-risk patients had elevated mortality and reduced survival time in comparison with their low-risk counterparts (Fig. 3E and F). PCA showed good discriminative performance of the necroptosis signature in TCGA and validation cohorts (Fig. 3G and H). Kaplan–Meier survival analysis with a two-sided log-rank test in the TCGA and validation cohorts showed patients in the high-risk group had significantly shorter OS compared to the patients in the low-risk group (P < 0.05) (Fig. 3I). The prognostic value of the signature was further tested in the validation cohort whereby a similar significant difference in OS was observed between the low- and high-risk groups (P < 0.05) (Fig. 3J). The AUC for the 3-year (0.712) and 5-year (0.711) survival rates in the TCGA cohort and the 3-year (0.703) and 5-year (0.740) survival rates in the validation cohort showed favorable specificity and sensitivity of the signature in predicting OS (Fig. 3K and L).

Necroptosis-related signature constructed and validated in the TCGA and GEO cohorts, repsectively. (A) The 7 necroptosis-related genes with P < 0.05 and their hazard ratios from univariate Cox proportional hazards regression analysis. (B) Tuning parameter selection in the LASSO model for OS-relevant genes in the TCGA cohort. (C) The LASSO coefficient profile of the 6 necroptosis-related genes in the TCGA cohort. The vertical line indicates the coefficient selected by LASSO. (D) Forrest plot showed that a total of 4 necroptosis-related genes were identified as prognosis-related by multivariate cox analysis. (E and F) The distribution and value of the risk scores in the TCGA and validation cohorts. (G and H) PCA of melanoma patients according to the risk score in the TCGA and validation cohorts. (I and J) Kaplan–Meier curves show that the low-risk group has significantly longer overall survival compared to the high-risk group in both the TCGA and validation cohorts. (K and L) ROC curve analysis shows the predictive efficiency of the established risk score in TCGA and validation cohorts.

Correlation between risk score and TME

We performed a CIBERSORT algorithm to evaluate the infiltrating level of immune cells in the tumor microenvironment and made comprehensive comparisons with the risk score. The high-risk score was negatively associated with infiltration of immune cells, including plasma cells, activated memory CD4 + T cells, CD8 + T cells, activated NK cells, M1 macrophages, regulatory T cells (Tregs), and memory B cells, and positively associated with infiltration of M0 macrophages, M2 macrophages, and neutrophils (Fig. 4A). Furthermore, we observed that the risk score had a strong negative correlation with the TME scores (Fig. 4B).

Immune landscape between low- and high-risk groups. (A) The infiltration levels of 22 immune cell types in the two risk groups. (B) The relationship between signature and the TME score.

Correlation between risk score and response to immunotherapy

We explored the relationships between risk score and response to immunotherapy and evaluated whether risk score could serve as a predictor to immune response. We selected 47 immune checkpoints to analyze their expression differences between the two risk subgroups and observed a higher expression of 41 immune checkpoints in patients with low-risk scores (Fig. 5A). In The Cancer Immunome Atlas, the IPS which was based on immunogenicity could achieve high accuracy in predicting the immunotherapy response of patients. Therefore, we analyzed the relationship of IPS between high and low-risk score groups. Interestingly, we found that the total IPS for PD-1 or CTLA-4, or both PD-1 or CTLA-4 blockers in the low-risk score group was significantly higher than that in the high score group, which meant that patients with lower risk scores tended to have a better response to immunotherapy (Fig. 5B–D). Next, the TIDEs were also calculated to predict response to immunotherapy. Similarly, patients with the low-risk score have a lower TIDE score (P < 0.05), suggesting a better response to immunotherapy (Fig. 5E). Taken together, the necroptosis signature established by us could serve as a robust predictor to immune response.

The relationships between risk score and immunotherapy. (A) Difference of immune checkpoint expression between high- and low-risk groups. (B–D) Difference of immunogenicity between high- and low-risk groups. (E) Difference of TIDE score between low- and high-risk groups. (F) Comparison of pigmentation score in the high- and low-risk groups. (G) Survival analysis of high and low pigmentation score group. (H) Survival analysis of melanoma stratified by pigmentation score and risk score.

Correlation of risk score with pigmentation score

We compared the pigmentation score of the high- and low-risk groups. The result indicated that the pigmentation score of the high-risk group was visibly higher than those of the low-risk group (Fig. 5F). Based on the median value of the patient’s pigmentation score, we divided the patients into low and high pigmentation score subgroups. KM survival analysis indicated that the OS rate in the high pigmentation score group was obviously lower than that in the low pigmentation score group (Fig. 5G). Stratified survival analysis demonstrated that the prognostic value of the risk score was not affected by the pigmentation score. Risk score subtypes showed significant differences in survival between low and high pigmentation score subgroups (Fig. 5H). Taken together, the risk score may be an effective indicator that is independent of the pigmentation score and can effectively evaluate the prognosis of patients.

Independent prognostic significance of signature and stratification analysis

Univariate analysis indicated that among the pertinent clinicopathological parameters, high age, advanced stage, and high-risk score were poor prognostic factors (Fig. 6A). In multivariable analysis, after adjustment for other confounders, age, stage, and the risk score were still independent prognostic indicators (Fig. 6B). To confirm the prognostic discriminatory power of the signature, we performed stratified survival analysis in various clinical subgroups, including age (age < 60 and age > 60), gender (female and male), tumor location (primary tumor, regional lymph node, and distant metastasis), stage (I-II and III-IV), Breslow depth (< = 1.5 and > 1.5), ulceration (no and yes), and tumor status (tumor-free and with tumor). As the result shown in Fig. 6C, the OS of the low-risk patients based on age (P < 0.001 in both age < 60 and age > 60), stage (P < 0.001 in both I-II and III-IV), tumor location (P = 0.005 in primary tumor; P < 0.001 in regional lymph node), Breslow depth (P = 0.005 in <  = 1.5, P < 0.001 in > 1.5), ulceration (P = 0.017 in no and P < 0.001 is yes), and tumor status (P < 0.001 in tumor-free and P = 0.007 in with tumor) was significantly higher than that of high-risk patients.

Univariate and multivariate Cox proportional hazards regression analysis for pertinent clinicopathological parameters. (A) Higher age, advanced stage, and high-risk groups were significantly poor prognostic factors in univariate analysis. (B) The high-risk group remains as an independent poor prognostic factor (P < 0.05) in the multivariate analysis. (C) Survival rates of the high- and low-risk patients in the subgroups based on clinicopathological characteristics.

Drug sensitivity analyses

We compared the sensitivity of two risk groups to common anticancer drugs to identify potential cutaneous melanoma treatment modalities. As shown in Supplementary Fig. 2, patients in the high-risk group were more sensitive to Docetaxel, Elesclomol, and Pazopanib, while those in the low-risk group were more sensitive to Gemcitabine, Gefitinib, Cyclopamine, Lenalidomide, and Metformin.

Construction of prognostic nomogram

To explore the potential value of the signature in clinical practice, we constructed a nomogram based on risk scores and clinical variables to predict the 3- and 5-year survival rates. A nomogram was constructed integrating the age, TNM stage, and risk score in the multivariate analysis (Fig. 7A). AUCs of time-dependent ROC curve at 3 and 5 year-OS for the nomogram-based prediction model were 0.726 and 0.704 (Fig. 7B). Calibration plots revealed that the nomogram showed perfect concordance between the observed and predicted survival rates at 3- and 5-year (Fig. 7C).

Construction and validation of the nomogram. (A) Nomogram for predicting 3-year and 5-year OS of melanoma patients incorporating independent prognostic factors from multivariate Cox proportional hazards regression. (B) AUCs of time-dependent ROC curves at 3- and 5-year OS for the independent prognostic factors. (C) Calibration curves corrected for deviations in agreement between the predicted and observed OS rates at 3 and 5 years.

Gene Set Enrichment Analysis (GSEA)

We conducted GSEA between the low- and high-risk subgroups using the entire gene network based on the TCGA dataset. We observed that several pathways with enrichment in the low-risk group were related to immunity with the filter criteria of FDR q value < 0.05 (Table S5), including chemokine signaling pathway, T cell receptor signaling pathway, Toll-like receptor signaling pathway, natural killer cell-mediated cytotoxicity, and B cell receptor signaling pathway (Supplementary Fig. 3).