Abstract Necroptosis is a novel programmed cell death that affects the tumor heterogeneity, microenvironment, and prognosis, which is not well elucidated in skin cutaneous melanoma (SKCM). The SKCM-TCGA, [32]GSE65904, and [33]GSE215120 datasets were downloaded from the TCGA and GEO databases, respectively. The necroptosis-related genes were identified by weighted co-expression network analysis (WGCNA) and single-cell sequencing analysis. COX and LASSO regression was used to construct the prognostic model. Survival analysis, immune infiltration analysis, and tumor mutation analysis between the high-necroptosis score (NCPTS) and low-NCPTS groups were performed. Finally, real-time PCR experiment was carried out to verify the results. A prognostic model based on 9 necroptosis-related genes (TUFM, CD53, CLEC2D, KLRC1, STAT4, IFI35, XCL1, TAPBP, and SOD2) was constructed to predict the prognosis of SKCM patients. The patients in high-NCPTS group had a poor prognosis. The expression of immune checkpoint-related gene and drug sensitivity were higher than those in the low-NCPTS groups, indicating susceptible for immunotherapy. Real-time PCR showed that TUFM expression was significantly higher in A375 cells than control (P < 0.05). Besides, TUFM expression was also validated by TISCH and HPA database. The prognostic model might provide guidance for the prognosis and immunotherapy for SKCM patients, contributing to a better understanding of necroptosis in SKCM. Supplementary Information The online version contains supplementary material available at 10.1038/s41598-025-07829-2. Keywords: Skin cutaneous melanoma, Necroptosis, Single-cell sequencing, Prognostic model, Immunotherapy Subject terms: Cancer, Cancer models, Skin cancer, Tumour immunology Introduction Skin cutaneous melanoma (SKCM), originating by uncontrolled growth of melanocytes, leads to most of the skin cancer deaths though it accounts for less than 5% of all skin cancers^[34]1. The vast majority of early-diagnosed cases are treated with curative-intent surgical resection, and the prognosis is favorable^[35]2. However, patients with resected stage II and III melanoma still have a greater risk of recurrence. For the recurrent, metastatic, and high-risk (thicker primary melanomas, ulceration, or more extensive lymph node involvement) patients, surgery is no longer sufficient and effective^[36]3. Therefore, chemotherapy, targeted therapy, adjuvant therapy, immunotherapy, and/or their combinations were of great importance for treatment^[37]4. Despite the major advances that have been achieved in improving outcomes for SKCM patients, the 5-year survival rate is less than 20% for metastatic patients^[38]1. Therefore, developing novel prognostic molecular markers and therapeutic agents or methods was urgent for SKCM patients’ survival. Necrosis is believed to be an “accidental” type of death, which is not regulated by molecular events and the opposite modality of cell death compared to apoptosis^[39]5. In 2015, necroptosis was first proposed as a novel programmed form of cell death, which harbored unique molecular mechanisms and was different from apoptosis and necrosis^[40]6,[41]7. Necroptosis is mainly mediated by RIPK1 (receptor-interacting protein kinase 1), RIPK3, and MLKL (mixed lineage kinase domain-like pseudokinase) and is characterized to be inhibited by the necrostatin-1 (Nec-1)^[42]6,[43]7. Necroptosis is of great significance for cancer because drug resistance to apoptosis is one of the major hallmarks of tumors. Nevertheless, necroptosis has been reported to be a friend or a foe of cancer. Necroptosis’ dual effects of promoting or reducing tumor growth have been reported in different types of cancer, such as colorectal cancer (CRC), hepatocellular carcinoma (HCC), or pancreatic ductal adenocarcinoma (PDAC)^[44]6–[45]8. In addition, several chemotherapy or immunotherapy drugs (including neoalbaconol, sorafenib, and obatoclax) can cause necroptosis to enhance antitumor immunity^[46]6,[47]7. Therefore, its specific role in cancer remains unclear, and the relationship between necroptosis and SKCM is worth exploring. In the present study, we conducted comprehensive bioinformatic analysis on TCGA and GEO datasets of SKCM patients, including single-cell sequencing analysis, expression analysis, survival analysis, and mutation analysis. Upon these analyses, we developed a necroptosis-related prognostic signature for SKCM patients. Then the SKCM patients were divided into high-necroptosis score (NCPTS) and low-NCPTS groups based on the median signature scores. The SKCM patients with high-NCPTS group showed a poor prognosis. In addition, tumor mutation load, immune microenvironment, and survival analysis were also explored. The immune infiltration levels and immune checkpoint-related genes were highly expressed in the high-NCPTS group compared to those in the low-NCPTS group. Finally, we verified the function of the TUFM gene using other databases. Overall, our study might provide new insights into the understanding of immune mechanism, diagnosis, and treatment of SKCM. Materials and methods Transcriptome data collection The clinical data, transcriptome data, and mutation data of SKCM patients were downloaded from The Cancer Genome Atlas (TCGA) database. After removing duplicates, non-tumor samples, and missing values, 452 TCGA-SKCM patients were finally included. The [48]GSE65904 dataset (array data) from the Gene Expression Omnibus (GEO) database was collected. The primary data were standardized and analyzed using R software (version 4.2.1) with limma package (version 3.52.4). Single cell sequencing data analysis The single-cell RNA sequencing (scRNA-seq) dataset [49]GSE215120 of SKCM, containing 11 samples in total, was downloaded from the GEO database. The [50]GSM6622295 and [51]GSM6622302 were excluded because of treatment and individual lymphocyte populations filter. The load and analysis of scRNA-seq data was achieved using R software with Seurat package (version 4.4.0). The criteria for cells and genes included in the study were as follows: (i) genes were expressed in at least 3 cells; (ii) cells with expressed genes fluctuated between 200 and 7000; (iii) cells with less than 10% mitochondrial genes were retained; (iv) doublets were removed. The number of highly variable genes was set at 3000. The “SCT” standardized method was used to integrate samples and remove the batch effect. The dimension was set as 20 and the UMAP algorithm was used to reduce the dimension and visualization of data. The k-nearest neighbor (KNN) algorithm was used for dimension reduction and clustering analysis. Next, the cells were annotated using the signatures from the original publication at a resolution of 2^[52]9,[53]10. The distribution and expression level of necroptosis-related genes in different cell populations were displayed subsequently. Sources of necroptosis-related genes and single sample gene set enrichment analysis 909 necroptosis-related genes were acquired from the GeneCards database ([54]https://www.genecards.org/), and 397 genes with a correlation greater than 0.5 were retained. Then the enriched scores associated with necroptosis in each SKCM patient were analyzed by single sample gene set enrichment analysis (ssGSEA) algorithm. Weighted gene co-expression network analysis Weighted gene co-expression network analysis (WGCNA) was used for detecting highly covarying gene sets, and candidate biomarker genes, as well as to investigate the associations between hub genes and immune phenotype. In the present study, the necroptosis-related gene set in SKCM was identified using WGCNA analysis. Clustering was performed using a soft threshold power of 10, a deepSplit of 2, and a minimum module size of 80. Construction and evaluation of the prognostic model Univariate Cox analysis was firstly carried out to screen necroptosis genes associated with prognosis. Then, the pivotal genes related to patient prognosis were identified using the least absolute and selection operator (LASSO) regression. Finally, a prognostic model was developed. All patients were divided into the high-NCPTS and low-NCPTS groups based on the median NCPTS calculated from the formula. Survival analysis was performed for the two groups and the accuracy of the model was also assessed. External validation of the model [55]GSE65904 was used as an external validation cohort to further evaluate the model. Similarly, patients in [56]GSE65904 cohort were divided into the high-NCPTS and low-NCPTS groups based on the median NCPTS and survival analysis was performed. Moreover, the receiver operating characteristic (ROC) curve and principal component analysis (PCA) were also performed to evaluate the accuracy of the model. Immune microenvironment and tumor mutation load analysis Seven different algorithms, including CIBERSORT, EPIC, Estimate, MCP-counter, Quantiseq, TIMER, and xCell, were used to assess the relationship between the risk score and immune filtration status for the patients in the TCGA database. CIBERSORT method are based on support vector regression using a linear kernel. MCP-counter is applied gene set enrichment approaches. Quantiseq and TIMER are based on linear approaches. ESTIMATE and xCell are based on log-linear regression for gene set enrichment analysis (GSEA). The two-way Pearson correlation was applied to compare the immune score between the high-NCPTS and low-NCPTS groups. The immune microenvironment (TME) analysis was using R software with the IOBR package (version 0.99.8). Besides, we also compared the differentially expressed immune checkpoint-related genes between the high-NCPTS and low-NCPTS groups. In addition, the SKCM mutation data were downloaded from the cBioPortal database ([57]https://www.cbioportal.org/datasets). The first 20 mutated genes in the high-NCPTS and low-NCPTS groups were presented, respectively. the tumor mutation load (TMB) of each sample was calculated using the function TMB of the maftools package (version 2.12.0) in R software. Function and pathway enrichment analysis The differently expressed genes (DEGs) between the high-risk and low-risk groups were identified with a threshold of P < 0.05 and log|FC (fold change)|>1. Gene Oncology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG)^[58]11,[59]12 analyses were performed to explore the potential signaling pathways of these DEGs using R software with the clusterProfiler package (version 4.4.4). OncoPredict for drug sensitivity analysis The OncoPredict R package was used to predict drug responses for cancer patients in vivo^[60]13. OncoPredict matches the gene expression profile of tissues to the half-maximal inhibitory concentration (IC50) of various cancer cell lines to drugs from the Genomics of Drug Sensitivity in Cancer (GDSC; [61]https://www.cancerrxgene.org/). A total of 198 drugs were calculated, and the drug sensitivities between the high-NCPTS and low-NCPTS groups were compared using Wilcoxon tests. A P < 0.05 was considered to be significant. Construction of the nomogram A nomogram was constructed based on NCPTS and clinical indicators to provide a quantitative tool to assess the individual survival risk of SKCM, which integrated the 1-, 3-, and 5-year overall survival (OS) of SKCM patients. The Kaplan-Meier survival analysis was carried out using the R software with the survminer package (version 0.4.9). The 1-, 3-, and 5-year OS predicting outcomes were using the R software with the survivalROC package (version 1.0.3). To evaluate the nomogram’s predictive performance, the calibration curve, concordance index (C-index), and decision curve were also performed. The verification of TUFM gene we explored TUFM gene expression in immune and non-immune cells in 10 SKCM single-cell datasets using the tumor immune single-cell hub (TISCH) database ([62]http://tisch.compgenomics.org). Then we used the Human Protein Atlas (HPA) database ([63]https://www.proteinatlas.org/) to verify whether TUFM expression in SKCM was different from that in normal skin tissues at the protein level. Finally, quantitative real-time polymerase chain reaction (qRT-PCR) was used to identify the TUFM expression in A375 cells and human PIG1 cell. A375 cells were cultured in DMEM (Mediatech lnc., VA, USA) supplemented with 10% fetal bovine serum (FBS, PAN, Seratech, Germany) and 1% penicillin-streptomycin (P/S) solution. The cells were placed in the humidified chamber of 37 °C with 5% CO[2]. PIG1 cells were cultured in RPMI-1640 (Gibco, Thermofisher lnc., USA) supplemented with 10% FBS (PAN, Seratech, Germany) and 1% P/S solution, which was also placed at 37℃ with 5% CO[2]. Total RNA Kit II (Code No. R6934-01, OMEGA, USA) was used to extract the total RNA of cell lines according to the manufacturer’s instructions. The RevertAid First Strand cDNA Synthesis Kit (Code No. K1622, Thermo Fisher Scientific, Waltham, MA, USA) was used for cDNA synthesis. qRT-PCR was performed using the Taq Pro Universal SYBR qPCR Master Mix (Code No. Q712-02, Vazyme, Nanjing, China) and amplified by cobas z480 (Roche Molecular Systems, Inc, New Jersey, USA). GAPDH was used as an internal control. All primers were synthesized by Sangon Biotech (Shanghai, China). The relative gene expression level was determined through the 2−∆∆Ct method. The qRT-PCR assay was repeated three times for each gene. The used primer sequences were as follows: TUFM (products:139 bp) forward 5′-TGTGGTGGTGTATGTGAAC-3′, reverse 5′-GAGAGCAGAGCCTACGAT-3′; GAPDH (products:101 bp) forward 5′-TCTCTGCTCCTCCTGTTC-3′, reverse 5′-GTTGACTCCGACCTTCAC-3′. Results Single cell sequencing data analysis 9 samples extracted from the [64]GSE215120 scRNA-seq dataset were integrated without an obvious batch effect and further analyzed (Fig. [65]1A). Finally, we obtained 26,287 cells containing 33,538 features after strict quality control and removing doublets. Then we used the UMAP algorithm to reduce the dimension and clustered all the cells into 31 clusters by KNN method (Fig. [66]1B and C). 397 necroptosis-related genes from the GeneCards database were input and the necroptosis score in each cell was calculated by using the PercentageFeatureSet function (Supplementary Table [67]S1). The cells were divided into the high-NCPTS and low-NCPTS groups compared to their median necroptosis score and were presented by the t-SNE algorithm (Fig. [68]1D). Next, eight distinct cell types of B cells, CD4^+T cells, CD8^+T cells, endothelial cells, fibroblasts, melanoma cells, monocytes and macrophages, NK cells, and T cells were annotated according to the surface marker genes of different cell types and their expression in different clusters (Fig. [69]1E and Supplementary Table [70]S2). Fig. 1. [71]Fig. 1 [72]Open in a new tab The flowchart of the study. Analysis of cell developmental trajectory and cell-cell communication networks for SKCM We used the monocle to construct the developmental trajectory of immune cells and melanoma cells (Fig. [73]2A). Results showed that the immune and melanoma cells exhibited various states (Fig. [74]2A). Figure [75]2B showed the over-expressed ligands or receptors and their interactions in all cell types to identify their interactions. Circle diagrams showed the integrated cell communication networks of interacting times and strengths between arbitrary two cell groups (Fig. [76]2C and D). Monocytes and macrophages were demonstrated to contribute the most in immune communication networks compared with other immune cells (Fig. [77]2C and D). Besides, Fig. [78]2E and F revealed the distinct contributive signals among different immune cell subgroups. Fig. 2. [79]Fig. 2 [80]Open in a new tab Single-cell sequencing analysis of the [81]GSE215120 dataset. (A) The integration effect of 9 samples was pretty good. (B and C) Dimensionality reduction and cluster analysis. All the cells in 9 samples were clustered into 31 clusters. (D) The cells were divided into the high-NCPTS and low-NCPTS groups. (E) The cells were annotated as B cells, CD4^+T cells, CD8^+T cells, monocytes and macrophages, NK cells, T cells, endothelial cells, fibroblasts, and melanoma cells, respectively. WGCNA analysis As mentioned, the covarying gene expression associated with necroptosis was performed using the WGCNA analysis for the TCGA cohort. We selected 10 as the soft threshold power with a R^2 > 0.8 and stable mean connectivity trends (Fig. [82]3A and B). A total of 16 non-grey modules were acquired with a minimum modules size of 80 and a deepSplit of 2. The deepSplit parameter is a key option in hierarchical clustering tree cutting, used to control the sensitivity of module recognition and affect the degree of module subdivision. The deepSplit of 2 represents the moderate granularity in the splitting of the dendrogram. Finally, a total of 9 non-gray modules were acquired with the similarity domain value of 0.4 (Fig. [83]3C). Then we selected the four non-gray modules MEgreenyellow, MEgreen, and MEsalmon for subsequent analysis (1963 genes), which were most closely associated with necroptosis (all P value < 0.0001, Fig. [84]3D). Fig. 3. [85]Fig. 3 [86]Open in a new tab Weighted gene co-expression network analysis (WGCNA). (A and B) WGCNA. When the soft domain value was 10, R^2 > 0.8, the data was consistent with the power-law distribution. The mean connectivity was stable and suitable for subsequent analysis. (C) A total of 16 non-grey modules were acquired with a minimum module size of 80 and a deepSplit of 2. Next, the similarity domain value was set as 0.4, and the modules lower than 0.4 were merged to acquire 9 non-gray modules. (D) The correlation between MEgreenyellow module and necroptosis was strongest (correlation coefficient = 0.82, P < 0.0001). Construction and validation of necroptosis-related prognostic model in the TCGA patients 1963 genes closely related to necroptosis were acquired for subsequent analysis by matching TCGA and [87]GSE215120 coexisting genes (Supplementary Table [88]S3). In the TCGA cohort, 194 genes related to the prognosis of SKCM patients were obtained by univariate COX analysis with P < 0.05. Then, 9 genes were included in the model by using LASSO regression analysis (Fig. [89]4A and B). The genes included in the model by LASSO regression analysis were summarized in Table [90]1. Finally, 9 genes were used to construct the prognostic model. The prognostic model was calculated as follows: NCPTS = TUFM*0.23039 + CD53*(-0.18404) + CLEC2D*0.10385 + KLRC1*(-0.461 49) + STAT4*(-0.01027) + IFI35*(-0.53473) + XCL1*(-1.16389) + TAPBP*0.0 0527 + CLEC2D*(-0.03740). Therefore, the patients were divided into the high-NCPTS and low-NCPTS groups based on the median value. As shown in Fig. [91]4C and D, patients with high-NCPTS value had poor prognosis compared to the low-NCPTS group in the TCGA and [92]GSE215120 cohorts, respectively. Then we conducted the ROC curve analysis to assess the accuracy of NCPTS for the prognosis of SKCM patients. Results showed that the AUCs of 1-, 2-, 3-, and 5-year were 0.622, 0.662, 0.653, and 0.665 in the TCGA cohort, respectively (Fig. [93]4E). In the [94]GSE215120 cohort, the AUCs of 1-, 2-, 3-, and 5-year were 0.639, 0.629, 0.651, and 0.636, respectively (Fig. [95]4F). Moreover, PCA analysis was performed on the model in the training set and validation set, respectively. we found that the model could distinguish the patients well in both the training cohort and the validation cohort (Fig. [96]4G and H). Fig. 4. [97]Fig. 4 [98]Open in a new tab Construction and validation of necroptosis-related prognostic model. (A and B) Nine genes were selected to construct the prognostic model using LASSO regression analysis. (C and D) The patients with high-NCPTS group had worse prognosis in the TCGA cohort (P < 0.0001) and [99]GSE65904 (P = 0.00051), respectively. (E) The AUCs of 1-, 2-, 3-, and 5-year for the model in the TCGA cohort were 0.622, 0.662, 0.653, and 0.665, respectively. (F) The AUCs of 1-, 2-, 3-, and 5-year for the model in the [100]GSE65904 cohort were 0.639, 0.629, 0.651, and 0.636, respectively. (G and H) PCA analysis showed that the model could group SKCM patients well in the TCGA cohort and [101]GSE65904 cohort, respectively. Table 1. Coefficients in the least absolute and selection operator (LASSO) model. ID Gene Coefficient 1 TUFM 0.23039 2 CD53 − 0.18404 3 CLE2D 0.10385 4 KLRC1 − 0.46149 5 STAT4 − 0.01027 6 IFI35 − 0.53473 7 XCL1 − 1.16389 8 TAPBP 0.00527 9 SOD2 − 0.03740 [102]Open in a new tab Immune infiltration analysis and mutation landscape The TME plays a vital role in the prognosis of patients. Immune infiltration analysis between the high-NCPTS and low-NCPTS groups was also evaluated. To provide a reference for immunotherapy, we explored the differences in immune infiltration levels among different groups. As shown in Fig. [103]5A, the Stromal score, Immune score, and Estimate score were higher in the high-NCPTS group, while the tumor purity was lower compared with the low-NCPTS group. Besides, CD4^+T cells, CD8^+T cells, Th cells, and macrophages showed more infiltrations in the high-NCPTS group than those in the low-NCPTS group (Fig. [104]5A). In addition, we also explored the expressions of genes related to immune checkpoints. Except for KDR and TGFBR1, all the immune checkpoint-related genes were highly expressed in the high-NCPTS group than those in the low-NCPTS group (Fig. [105]5B). The tumor cells could escape anti-tumor immunity through an immunosuppressive mechanism of immune checkpoints. Therefore, the patients of the high-NCPTS group might have more benefits in treating them with immune checkpoint inhibitors (ICI). Next, we displayed the top 20 mutated genes in the high-NCPTS and low-NCPTS groups using mutation landscape analysis. TTN, MUC16, BRAF, and DNAH5 were the 4 top mutated genes in both the two groups (Fig. [106]5C and D). Fig. 5. [107]Fig. 5 [108]Open in a new tab Immune infiltration analysis and mutation landscape. (A) Heatmap of immune cell infiltration in the high-NCPTS and low-NCPTS groups. (B) Comparison of immune checkpoint-related genes between the high-NCPTS and low-NCPTS groups. (C and D) The mutation landscape between the high-NCPTS and low-NCPTS groups in the TCGA cohort. **P < 0.01; ****P < 0.0001; ns, no significance. Enrichment analysis and drug sensitivity Figure [109]6A showed that DEGs were mostly enriched in immune-related pathways indicated by GO analysis, such as leukocyte mediated immunity, positive regulation of cell activation, positive regulation of leukocyte activation, and immune response-regulating signaling pathway (Supplementary Table [110]S4). KEGG analysis revealed that DEGs were mostly enriched in the cytokine-cytokine receptor interactions, chemokine signaling pathway, and cell adhesion molecules (CAMs) (Fig. [111]6B and Supplementary Table [112]S4). Besides, we compared the drug sensitivity in the two risk groups based on the GDSC screening data. The IC50s of cisplatin, cyclophosphamide, dabrafenib, and trametinib were significantly different (Fig. [113]6C–F). SKCM patients in the high-NCPTS group were more sensitive to those chemotherapy and targeted drugs (Fig. [114]6C–F). Fig. 6. [115]Fig. 6 [116]Open in a new tab Function and pathway enrichment analysis and drug sensitivity analysis. (A) GO analysis of DEGs. (B) KEGG analysis of DEGs. (C–F) The comparisons of the IC50s of cisplatin, cyclophosphamide, dabrafenib, and trametinib between the high-NCPTS and low-NCPTS groups. Cell localization of 9 model genes scRNA-seq was performed to evaluate the cell localization of model genes. As shown in Fig. [117]7A and B, TUFM, IFI35, and TAPBP were mainly expressed in tumor cells, CD53 and TAPBP were mainly expressed in B cells, CLEC2D was mainly expressed in CD4^+T cells, SOD2 was mainly expressed in monocyte and macrophages. TUFM and TAPBP were expressed in all cell types. Fig. 7. [118]Fig. 7 [119]Open in a new tab The cell localization of the 9 model genes by single-cell sequencing analysis. (A and B) TUFM, IFI35, and TAPBP were mainly expressed in tumor cells, CD53 and TAPBP were mainly expressed in B cells, CLEC2D was mainly expressed in CD4^+T cells, SOD2 was mainly expressed in monocyte and macrophages. Construction of the nomogram Intending to better evaluate the prognostic prediction for SKCM patients, we constructed a nomogram containing NCPTS and clinical features, including age, gender, T stage, and stage. The predicted 1-, 3-, and 5-mortality rates of SKCM patients was 5.15%, 26.3%, and 39.2%, respectively (Fig. [120]8A). Besides, the calibration curve indicated the consistency between the predicted OS and actual observed OS (Fig. [121]8B), which might be applied in clinical settings and guide making decisions. The decision curve analysis showed that the 5-year net benefit was larger than those of 1-year and 3-year, suggesting the nomogram a good effect in predicting the prognosis of SKCM patients (Fig. [122]8C). In addition, the C-index of the nomogram was higher than other clinical indicators (Fig. [123]8D). suggesting the prognostic ROC was also performed. Figure [124]8E showed the 1-, 3-, and 5-year AUC of the nomogram prediction were 0.64, 0.70, and 0.71, respectively, indicating good accuracy in predicting the survival of SKCM patients. Fig. 8. [125]Fig. 8 [126]Open in a new tab Construction of the nomogram model. (A) Nomogram predicting the 1-, 3-, and 5-year mortality rates for SKCM patients were 5.15%, 26.3%, and 39.2%, respectively. (B) Calibration curves of the nomogram. (C) The net benefit rate of 1-, 3-, and 5-year predicted the prognosis of SKCM patients. (D) The time-dependent C-index curves. (E) 1-, 3-, and 5-year ROC curves predicting the performance of the nomogram. Survival analysis and validation of TUFM The SKCM patients with high expression of TUFM showed a worse prognosis than those with low expression of TUFM (P < 0.0001, Fig. [127]9A). In addition to TUFM, the high expression of the other eight genes also indicates a worse prognosis for SCKM patients (Supplementary Fig. 1). Besides, the TUFM expression was significantly higher in A375 cells compared to the PIG1 cell (P < 0.001, Fig. [128]9B). Next, we evaluated the TUFM expression situations both in immune and non-immune cells in ten SKCM single-cell datasets based on the TISCH database. Results showed that TUFM was highly expressed in CD4^+T cells, CD8^+T cells, dendritic cells (DC), and melanoma cells in the microenvironment (Fig. [129]9C). Moreover, we found that the protein expression of IFITM3 in SKCM tissue was also higher than that in normal skin tissue using the immunohistochemical data of the HPA database (Fig. [130]9D). Fig. 9. [131]Fig. 9 [132]Open in a new tab Validation of TUFM expression in melanoma. (A) Survival analysis of high-expression and low-expression TUFM in TCGA-SKCM patients. (B) The TUFM expression in A375 cell and PIG1 by qRT-PCR assay. (C) Locating the expression situations of TUFM in immune and nonimmune cells in ten SKCM single-cell datasets from the TISCH database. (D) TUFM protein expression in normal and melanoma tissues using immunohistochemical staining. ***P < 0.001. Discussion A high rate of metastasis is an important characteristic of SKCM^[133]14. Although the treatment of metastatic SKCM has been revolutionized by the application of immune and targeted therapies in the last decades, especially single-agent programmed death 1 (PD-1) inhibitors and cytotoxic T-lymphocyte antigen 4 (CTLA-4) inhibitors. However, only a few of patients can benefit from immune checkpoint inhibitors (ICIs) therapy. In addition, the toxicity, expensive price, evaluation of medication indication, and unavailability have also compromised the prognosis. Therefore, a thorough understanding of the pathogenesis of SKCM and its interaction with the TME may provide new ideas and targets for the treatment of SKCM. With the increased understanding mechanism, the role of necroptosis in tumors has received much attention. Evidence has indicated that necroptosis plays a vital role in tumor initiation, progression, metastasis, and anti-cancer therapy by bypassing acquired or intrinsic apoptosis-resistance. Meanwhile, several compounds and multiple chemotherapeutic agents have been investigated in a variety of cancer cell lines and animal models by the mechanism of triggering necroptosis^[134]5. However, the role of necroptosis in SKCM and the concomitant therapy remains unclear. In the present study, we presented an extensive single-cell transcriptome analysis elucidating the necroptosis phenomenon in SKCM, which facilitated an in-depth comprehension of the TME and cell type associated with necroptosis. The scRNA-seq technology is rapidly developing, enabling the profiling of biological tissue at unprecedented genomic depth. The technology has allowed for comprehensive analysis of the immune system and cells in TME, though the composition and the infiltration degree of immune cells in different tumor types are different^[135]15. The interaction between immune cells and the interaction between immune cells and tumor cells plays important roles in the development, progression, and metastasis of tumors. Therefore, elucidating these communication processes between immune cells and tumor cells are vital for exploring effective anticancer immunotherapy strategies. The scRNA-seq measures the whole transcriptome at a single-cell level and distinguishes different cell types in tumor tissue, which overcomes the limitations of traditional RNA sequencing methods^[136]16. In our study, we displayed the distribution of necroptosis in the immune cell region, such as CD4^+T cells, CD8^+T cells, and monocytes and macrophages. Besides, monocytes and macrophages were found to be a dominant contributor to immune cell-cell communications in the TME of SKCM. Next, we identified necroptosis-related genes using WGCNA and scRNA-seq data. The prognostic signature of necroptosis-related genes was developed by univariate COX and LASSO regression. Therefore, the SKCM patients were divided into high-risk and low-risk groups based on the NCPTS. Obviously, the patients in the high-risk group had a worse prognosis. Currently, escaping from the surveillance of the immune system and metastasis are still the focus of tumor research. The application of immunotherapy for treating advanced SKCM has achieved preliminary advances in the past decades, such as PD-1 inhibitors and CTLA-4 inhibitors^[137]4,[138]17. However, developing other immune checkpoints for inhibiting T-cell activation and new combinations of checkpoint inhibitors were necessary for improving patients’ prognosis. In our study, the patients in high-NCPTS groups showed higher LAG3 expression compared to the low-NCPTS groups. Importantly, the LAG3 expression in high-NCPTS groups was also higher than the expression of PDCD1 (PD-1) and CTLA-4. Therefore, when the two inhibitors targeting PD-1 and CTLA-4 were less effective or ineffective, the drugs targeting LAG3 might be considered an alternative for treating advanced SKCM^[139]18. Besides, several studies have explored the combinations of immunotherapy with other treatments (such as radiotherapy, chemotherapy, and targeted therapy) to treat solid tumors, including SKCM^[140]4,[141]19,[142]20. As expected, the combination of immunotherapy and necroptosis-related therapy might be a potential direction. The present study showed the discrepancy of immune infiltration analysis, immune checkpoints, and tumor mutation load between the two groups, which might provide new insights into the combination of immunotherapy and necroptosis-based therapy for SKCM. In addition, we compared the drug sensitivity between high-NCPTS and low-NCPTS groups using the GDSC2 data. The SKCM patients with high-NCPTS were more susceptive to the chemotherapy and targeted drugs, suggesting a promising direction for SKCM. Several studies have previously reported the relationship between necroptosis and SKCM. Ma et al. conducted a pan-cancer analysis of necroptosis‑related gene signatures for predicting the prognosis of 33 tumors^[143]21. The necroptosis-related signature could predict clinical outcomes and distinguish immune infiltration status for SKCM^[144]21. Niu et al. developed a prognostic model containing 4 necroptosis-related genes (BOK, CD14, CYLD, and FASLG) to reveal the possible connections between necroptosis and melanin pigmentation using the TCGA and GEO databases^[145]22. Similarly, another study also constructed a prognostic model based on 10 necroptosis-related genes to evaluate the prognosis and immunotherapy response for SKCM using the TCGA and GEO databases^[146]23. Cao et al. found that FASLG, PLK1, EGFR, and TNFRSF21 among the necroptosis-related genes predicted the prognosis of SKCM patients well^[147]24. However, the necroptosis-related genes included in our model were different from the aforementioned studies. We selected 9 necroptosis-related genes (TUFM, CD53, CLEC2D, KLRC1, STAT4, IFI35, XCL1, TAPBP, and SOD2) to construct the prognostic nomogram. TUFM, a nuclear-encoded mitochondrial protein, was reported in several tumors, including lung cancer, esophageal cancer, gastric cancer, and so on^[148]25–[149]27. As shown in Fig. [150]7, TUFM expression was highest among the 9 model genes. Dysregulation of TUFM has been linked to tumor growth, metastasis, and therapeutic resistance^[151]25–[152]27. However, the role of TUFM in SKCM has not been discovered. We first identified the TUFM gene for its highest hazard ratio (HR) among the 9 model genes. The SKCM patients with higher TUFM expression had a worse prognosis. Besides, we also found that the TUFM expression was upregulated in SKCM cell line (A375) and SKCM tissues compared to normal control. Nevertheless, the present study had several limitations. First, the underlying roles of necroptosis in SKCM still need to be investigated. Second, the results should be further validated in future studies by a series of mechanistic and relevant animal experiments. In summary, our study constructed a prognostic model based on 9 necroptosis-related genes for SKCM by integrating single-cell sequencing, TCGA, and GEO datasets. We conducted an in-depth analysis and obtained the results of single-cell sequencing data, such as cell trajectory, cell communication, and signaling pathway enrichment. Most importantly, we validated the TUFM gene expression using cell lines and the HPA database. The signature could be used to predict the prognosis of SKCM patients and provide a reference for treating SKCM linked to necroptosis. Electronic supplementary material Below is the link to the electronic supplementary material. [153]Supplementary Material 1^ (54KB, csv) [154]Supplementary Material 2^ (26KB, xls) [155]Supplementary Material 3^ (8.6MB, xls) [156]Supplementary Material 4^ (599.5KB, xls) [157]Supplementary Material 5^ (17KB, docx) [158]Supplementary Material 6^ (195.7KB, docx) Author contributions Xin Fan: Formal analysis; Methodology; Project administration; Resources; Writing-original draft; Writing-review & editing. Ling-Yi Lu: Formal analysis; Funding acquisition; Software; Validation; Writing-review & editing. Si-Han Wang: Data curation; Validation; Visualization; Writing-review & editing. Qiong-Yan Zhou: Data curation; Investigation; Validation; Writing-review & editing. Bing-Jiang Lin: Conceptualization; Methodology; Resources; Supervision; Visualization; Writing-review & editing. All authors contributed to the study conception and design. All authors read and approved the final manuscript. Funding This work was supported by the Medical Health Science and Technology Project of Zhejiang Province (Grant NO. 2025KY1338 and 2022PY022). Data availability The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors upon reasonable cause. The single cell RNA-seq and array datasets used in the study were all downloaded from GEO database ([159]https://www.ncbi.nlm.nih.gov/geo/). The used TCGA dataset was downloaded from GDC Data Portal ([160]https://portal.gdc.cancer.gov/). The necroptosis-related genes were downloaded from GeneCards database ([161]https://www.genecards.org/). The immunohistochemical data was acquired from the HPA database ([162]https://www.proteinatlas.org/). Declarations Competing interests The authors declare no competing interests. Footnotes Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. References