Necroptosis-associated classification combined with tumor microenvironment characteristic analysis of cutaneous melanoma
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).
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).
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).
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.
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.
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).
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).