Clarifying new molecular subtyping and precise treatment of melanoma based on disulfidptosis-related lncRNA signature

Abstract

Disulfidptosis, the newest form of programmed cell death, is closely associated with the immune microenvironment of cancer cells. Long non-coding RNA (lncRNA) has also been found to play a crucial role in melanoma. However, the role of disulfidptosis-related lncRNA in melanoma remains unclear. Through bioinformatic analysis of the transcriptional, clinical, and pathological data from the TCGA-SKCM (The Cancer Genome Atlas—Skin cutaneous melanoma) database, we established a 2-Disulfidptosis-related lncRNA (DRL) prognostic model and a novel molecular subtype for melanoma. The survival and ROC curves of the 2-DRL prognostic model demonstrated its strong efficacy in predicting the prognosis of melanoma. The high-risk group of melanoma exhibited a significant decrease in ESTIMATEScore, ImmuneScore, and StromalScore, indicative of pronounced immune suppression and exhaustion. Subgroup C2 of melanoma displayed an immune-activated state, while subgroups C1 and C3 showed immune suppression and exhaustion, potentially leading to poorer prognosis. Subgroup C1 demonstrated better sensitivity to Zoledronate, UMI-77, Nilotinib, and Cytarabine. Subgroup C2 exhibited greater sensitivity to Ribociclib, XAV939, Topotecan, and Ruxolitinib. Subgroup C3 showed higher sensitivity to VX-11e, Ulixertinib, Trametinib, and Afatinib. This study revealed the immune microenvironment status and targeted drug sensitivity in melanoma patients with different risk scores and molecular subtypes, offering valuable guidance for clinical treatment and identifying significant DRL targets for future in-depth research.

Introduction

Cutaneous melanoma accounts for 75% of skin cancer-related deaths, and its global incidence has been steadily increasing [1]. Malignant melanoma is one of the most immunogenic tumors due to its incredibly high genomic mutation load, making it highly likely to elicit a specific adaptive anti-tumor immune response [2, 3]. The plasticity of melanoma cells leads to a phenomenon known as “immune escape,” in which cancer cells acquire a less immunogenic phenotype and possess the ability to suppress anti-tumor immune cells within the tumor microenvironment (TME) [4, 5]. Despite many factors contributing to immune escape being elucidated, the complete overhaul of therapeutic anti-tumor immune reinfusion strategies has not been achieved, resulting in persistently high mortality rates [6]. Clarifying how melanoma reprograms immune cells and the TME holds promise for innovative improvements in immune therapy, potentially rescuing anti-tumor immunity in melanoma.

Disulfidptosis is a newly discovered form of programmed cell death [7]. Under conditions of glucose starvation, cells with high levels of SLC7A11 experience reduced intracellular NADPH levels and cysteine accumulation, leading to aberrant binding of disulfide bonds to the actin cytoskeleton and ultimately causing cell collapse. The study also demonstrates that glucose transporter inhibitors can induce disulfide sagging in tumor cells, inhibiting tumor progression. Beyond SLC7A11, existing research has indicated that Disulfidptosis-related genes (DRGs) could potentially serve as therapeutic targets for melanoma [8, 9]. Targeting Disulfidptosis-related genes may offer an effective approach for treating melanoma.

Long non-coding RNAs (lncRNAs) are transcripts of over 200 nucleotides that are not translated into proteins. The encoding sites of lncRNAs are among the most numerous regulatory and functional units in the non-coding regions of the genome [10]. They play a crucial role in regulating gene expression and protein function by interacting with DNA, RNA, and proteins [10]. Research in other cancer types suggests that lncRNAs may regulate Disulfidptosis-Related Genes (DRGs), thereby affecting the immune microenvironment and tumor progression [11,12,13]. However, the specific role of Disulfidptosis-Related lncRNAs (DRLs) in melanoma is currently unclear. Further investigation may help to elucidate the function and mechanisms of DRLs in melanoma.

In this study, we utilized bioinformatics analysis of the TCGA-SKCM (The Cancer Genome Atlas-Skin cutaneous melanoma) database to investigate the transcriptome, clinical, and pathological information of melanoma. We established a two-DRLs prognostic model and a novel molecular classification of melanoma, revealing the immune microenvironment status and targeted drug sensitivity in melanoma patients with different risk scores and molecular subtypes. Our findings provide effective guidance for clinical treatment and valuable DRL targets for future in-depth research in melanoma.

Materials and methods

Transcriptomic data and clinical information collection from melanoma patients

We downloaded the transcriptome data of 59428 genes for 472 melanoma patients, along with paired non-cancerous tissue transcriptome data, from the official TCGA data portal at https://portal.gdc.cancer.gov/(8). These data were further analyzed based on the gene count numbers. Additionally, we obtained the tumor mutational burden (TMB) data for each patient, as well as complete clinical information for 470 patients, including futime, fustat, age, gender, stage, T stage, M stage, and N stage. The transcriptome data, TMB data, and clinical data were merged for further analysis.

Identification of differentially expressed DRLs in SKCM

Using R version 4.2.1 from https://cran.r-project.org/ and the latest version of RStudio from https://posit.co/download/rstudio-desktop/, we merged the transcriptome data, TMB data, and clinical data for further analysis. From previous disulfidptosis-related experimental articles (references: [8, 9, 13], we meticulously identified 14 genes (ACTB, CD2AP, ACTN4, DSTN, CAPZB, FLNB, FLNA, INF2, MYH10, IQGAP1, PDLIM1, MYL6, MYH9 and TLN1) associated with disulfidptosis (DRGs) which were clearly verified by previous experiments. We performed Pearson correlation analysis between 14 DRGs and all lncRNAs in TCGA-SKCM using the R limma package to identify differentially regulated lncRNAs (DRLs) (corFilter = 0.4 and pvalueFilter = 0.001 are considered DRLs) [14]. Additionally, we conducted differential analysis on the DRLs using the R limma package, comparing the non-cancerous tissues to cancer tissues (logFCfilter = 1), to identify DRLs that are differentially expressed in SKCM.

Construction and identification of 2-DRL prognostic models

Using the previously merged transcriptome and clinical data of SKCM, we performed univariate Cox regression analysis using the R survival package to identify differentially regulated lncRNAs (DRLs) associated with melanoma prognosis [13]. To avoid overfitting of the prognostic model, we conducted lasso regression analysis using the R glmnet package [14]. And we set maxit = 1000 (number of iterations) for cross-validation to ensure it works.

Based on the Cox model constructed using the R survival package, we identified a 2-DRL prognostic model with a risk score formula of (−1.31020221634025 * USP30-AS1) + (LINC02560 * 0.258515573577887). To visualize the expression levels of each gene in SKCM samples, we used the R pheatmap package to generate a heatmap [8]. For further analysis, we utilized GEPIA (http://gepia.cancer-pku.cn/) to analyze the expression levels and survival curves of USP30-AS1 and LINC02560 in SKCM. Subsequently, we validated the 2-DRL prognostic model by performing survival analysis and generating Kaplan–Meier survival curves using the R survival package and R survminer package. To assess the effectiveness of the model, we employed the R timeROC package to construct ROC curves.

Clarify the advantages of the 2-DRL prognostic model

Using R survival, we conducted univariate and multivariate Cox regression analyses to assess the superiority of the 2-DRL prognostic model in predicting SKCM prognosis compared to other clinical pathological factors. Through R regplot and R rms packages, significant clinical pathological factors identified via multivariate Cox regression were used to construct forest plots, evaluating the prognosis of melanoma patients at 1, 3, and 5 years [14]. Additionally, SKCM patients were stratified into various clinical subgroups based on clinical pathological factors to evaluate whether the 2-DRL model possesses similar prognostic diagnostic capabilities across different melanoma clinical subgroups [14].

Gene functional enrichment analysis

Using the R limma package, we performed differential analysis between the high-risk and low-risk groups, identifying significantly differentially expressed genes in both groups (logFC filter = 1, fdrFilter = 0.05). We use wilcoxTest to calculate the corresponding p-value, and then use the p.adjust function in R language (p.adjust (as.numeric (as.vector (pvalue)), method = "fdr")) to calculate the corresponding FDR value. To visualize the differential expression, we utilized the R pheatmap package to generate a heatmap and volcano plot. For the differentially expressed genes, we conducted enrichment analysis including Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Gene Set Enrichment Analysis (GSEA) using various R packages such as R clusterProfiler, R org.Hs.eg.db, R enrichplot, R ggplot2, R circlize, R RColorBrewer, R dplyr, R ggpubr, and R ComplexHeatmap [13].

Immune microenvironment analysis

The R limma and R estimate packages were employed to calculate the immune microenvironment scores and immune cell content in each melanoma sample. Specifically, the R limma and R ggpubr packages were utilized to compute the ESTIMATEScore, ImmuneScore, and StromalScore in melanoma samples. Additionally, the R limma, R survival, and R survminer packages were employed to analyze the relationship between immune cell content and prognosis in melanoma patients. To further investigate immune cell infiltration, immune function, and immune checkpoint gene analysis in melanoma samples, various R packages including R GSVA, R limma, R GSEABase, R ggpubr, and R reshape2 were utilized [13].

Drug sensitivity analysis

The R packages limma, oncoPredict, and parallel were utilized to calculate the IC50 (half maximal inhibitory concentration) values of 198 targeted drugs for each melanoma sample, serving as a measure of drug sensitivity [13]. Specifically, the lower the IC50, the more sensitive it is to drug treatment. We calculated the IC50 values ​​of each sample in different risk groups for different drugs to obtain the therapeutic drugs that are sensitive to patients in different risk groups.

TMB-related analysis

The R packages ggpubr and reshape2 were employed to calculate the TMB (Tumor Mutational Burden) in each melanoma sample and perform detailed analysis of specific mutated genes. By combining risk scoring with TMB, the R packages survival and survminer were used to more accurately predict the prognosis of melanoma patients [13].

Melanoma molecular classification

The R packages limma and ConsensusClusterPlus were utilized to perform unsupervised clustering analysis on melanoma samples, resulting in molecular subtyping. Cumulative Distribution Function (CDF) curves, triangle area plots, and clustering heatmaps were employed to determine the optimal molecular subtypes for melanoma samples. To validate the accuracy of the molecular subtyping, principal component analysis (PCA) was carried out using the R package ggplot2 [13].The analysis involved a comprehensive exploration of the molecular features of melanoma samples, aiding in better understanding the classification and characteristics of melanoma.

Results

Identification of lncRNAs associated with disulfidptosis in melanoma

We downloaded transcriptomic data of 59428 genes from 472 melanoma patients in the TCGA-SKCM dataset. Among these patients, 470 had complete clinical and pathological information. Subsequently, we extracted the transcriptomic data of all lncRNAs in SKCM, as well as 14 genes previously identified to play important roles in disulfidptosis. Using the R package limma, we performed Pearson correlation analysis between 14 DRGs and all lncRNAs in TCGA-SKCM, resulting in the identification of 83 lncRNAs (referred to as DRLs) significantly associated with disulfidptosis (corFilter = 0.4, pvalueFilter = 0.001) (Table S1). Next, we performed differential analysis on the identified 83 DRLs in SKCM using the R package limma, leading to the identification of 52 DRLs (logFCfilter = 1) exhibiting differential expression in SKCM (Table S2). For visualization, we used the R package pheatmap to generate a heatmap depicting the expression levels of the top 50 DRLs with the largest fold changes in SKCM samples (Fig. 1A). Additionally, we plotted a volcano plot to display the fold change of DRLs in SKCM (Fig. 1B).

Fig. 1

Revealing the identification of lncRNAs associated with prognosis and Disulfidptosis in melanoma. A Heatmap displays the expression of the top 50 Disulfidptosis-related lncRNAs in melanoma samples. B Volcano plot illustrates the differential expression fold change of Disulfidptosis-related lncRNAs in melanoma samples. C Heatmap represents the expression of differentially expressed Disulfidptosis-related lncRNAs in melanoma samples. D Forest plot demonstrates the association between Disulfidptosis-related lncRNAs and melanoma prognosis. E Regression coefficient pathway diagram displays the relationship between the number of included lncRNAs and the model fit in the study. F Cross-validation curve illustrates the relationship between the number of included lncRNAs and the model fit in the study. G Expression profiles of USP30-AS1 and LINC02560 in melanoma samples and paired adjacent non-cancerous samples based on the GEPIA website. H Survival curve shows the relationship between USP30-AS1, LINC02560, and melanoma prognosis

Using DRLs to build a prognostic model for melanoma

We endeavored to construct a prognostic model for melanoma using DRLs. By integrating transcriptomic and clinical data in SKCM, we conducted single-factor COX regression analysis using the R package survival, which resulted in the identification of 28 DRLs significantly associated with melanoma prognosis (coxPfilter = 0.05) (Table S3). Subsequently, we visualized the expression levels of these DRLs in SKCM samples using the R package pheatmap (Fig. 1C) and displayed the Hazard ratios of the DRLs in a forest plot (Fig. 1D). Next, we utilized the R package glmnet to perform lasso regression analysis to address overfitting in the prognostic model. The regression coefficient path diagram in Fig. 1E and the cross-validation curve in Fig. 1F indicated that the model achieved the best fit when incorporating two DRLs. By constructing a COX model, we determined a 2-DRL prognostic model (riskscore = (−1.31020221634025*USP30-AS1) + (LINC02560*0.258515573577887)). Furthermore, using the GEPIA website, we conducted expression and survival analysis of USP30-AS1 and LINC02560, revealing that USP30-AS1 is highly expressed in SKCM and associated with favorable prognosis (USP30-AS1 may be a tumor suppressor gene highly expressed in cancer tissues), whereas LINC02560 is highly expressed in SKCM and associated with poor prognosis (Fig. 1G–H).

Examination and validation of 2-DRL prognostic model

We validated the 2-DRL model and observed that high-risk and low-risk patients were evenly distributed across the All cohort, Train cohort, and Test cohort (Fig. 2A–C). Similarly, the distribution of deceased and surviving patients was also uniform across the All cohort, Train cohort, and Test cohort (Fig. 2D–F). Furthermore, the expression levels of USP30-AS1 and LINC02560 exhibited similar distributions across the All cohort, Train cohort, and Test cohort (Fig. 2G–I). Subsequently, we conducted survival analysis and constructed Kaplan–Meier survival curves using the R package survival and R package survminer. Remarkably, the high-risk score effectively predicted adverse prognosis for melanoma in the All cohort, Train cohort, and Test cohort (Fig. 2J–L). Additionally, ROC analysis using the R package timeROC indicated that the 2-DRL prognostic model achieved diagnostic efficiencies of 0.661, 0.664, and 0.672 at 1, 3, and 5 years respectively in the All cohort; 0.673, 0.664, and 0.684 at 1, 3, and 5 years, respectively, in the Train cohort; and 0.654, 0.658, and 0.656 at 1, 3, and 5 years, respectively, in the Test cohort.

Fig. 2

Testing and validation of the 2-DRL prognostic model for melanoma. A-C Line plots display the distribution of risk scores in melanoma samples from the all cohort (A), train cohort (B), and test cohort (C). D-F Dot plots illustrate the distribution of survival status in melanoma patients from the all cohort (D), train cohort (E), and test cohort (F). G-I Heatmaps demonstrate the expression profiles of USP30-AS1 and LINC02560 in melanoma patients from the all cohort (G), train cohort (H), and test cohort (I). J-L Survival curves show the prognostic status of melanoma patients with different risk scores in the all cohort (J), train cohort (K), and test cohort (L). M–O ROC curves demonstrate the effectiveness of the prognostic model in predicting the prognostic outcomes of melanoma patients in the all cohort (M), train cohort (N), and test cohort (O)

Next, we evaluated the advantages of the 2-DRL prognostic model compared to other clinical-pathological factors in predicting the prognosis of melanoma patients. We performed single-factor COX regression analysis using R survival to assess the hazard factors of age, stage, T stage, M stage, and riskscore in melanoma (Fig. 3A). The results showed that all these factors were significant hazard factors for melanoma, with hazard ratios greater than 1. In the multivariate COX regression analysis, age, T stage, M stage, and riskscore were identified as independent hazard factors for melanoma, and riskscore demonstrated the highest hazard ratio (1.923) among them (Fig. 3B). Furthermore, ROC curves revealed that riskscore had the highest diagnostic performance compared to other clinical-pathological factors (Fig. 3C). To further evaluate the prognostic ability, we employed R regplot and R rms packages to construct forest plots based on the clinically significant factors identified from the multivariate COX regression analysis. These plots assessed the 1-, 3-, and 5-year prognoses of melanoma patients (Fig. 3D). The forest plots indicated moderate effectiveness (C-index = 0.714). Additionally, we assessed the predictive capabilities of the 2-DRL model in various clinical subgroups of melanoma. The results demonstrated that the 2-DRL model significantly differentiated high-risk and low-risk groups in clinical subgroups such as Male, Female, < 65 years, > 65 years, stage 0–2, stage 3–4, T0-2, T3-4, N0-1, N2-3, and M0 (Fig. 3F–P). It is worth mentioning that although the high-risk group had worse prognosis in the M1 subgroup, it did not reach statistical significance, which could be attributed to the smaller sample size in the M1 clinical subgroup (Fig. 3Q). These analyses revealed the effectiveness and accuracy of the 2-DRL prognostic model in predicting the prognosis of melanoma patients in different clinical subgroups.

Fig. 3

The relationship between 2-DRL prognostic model and clinicopathological factors of melanoma. A Forest plot displays the results of univariate COX regression analysis of risk scores and other clinical-pathological factors in melanoma. B Forest plot displays the results of multivariate COX regression analysis of risk scores and other clinical-pathological factors in melanoma. C ROC curve demonstrates the predictive ability of risk scores and other clinical-pathological factors in predicting the prognosis of melanoma patients. D Line plot illustrates the construction of a melanoma diagnostic model based on clinically relevant factors with prognostic significance. E C-index curve reflects the accuracy of the constructed diagnostic model. F-Q Survival curves reflect the ability of the prognostic model to predict the prognosis of melanoma patients in different clinical subgroups

Melanoma is closely related to immune cells and immune function

We attempted to identify the distinct features of the high-risk and low-risk groups. To achieve this, we conducted differential analysis using the R limma package, which identified a total of 238 differentially expressed genes (logFC filter = 1, FDR filter = 0.05) (See Table S4). Subsequently, we utilized the R pheatmap package to visually represent the top 20 genes through a heatmap (Fig. 4A), and created a volcano plot to illustrate the fold changes of the differentially expressed genes (Fig. 4B). Further analysis of these differentially expressed genes using KEGG pathways revealed enrichment primarily in immune system-related functions, including Cytokine-cytokine receptor interaction, Cell adhesion molecules, Phagosome, Th1 and Th2 cell differentiation, and Th17 cell differentiation (Fig. 4C). Similarly, the GO analysis indicated enrichment in Biological processes such as activation of immune response, immune response-regulating signaling pathway, and immune response-activating signaling pathway. In terms of Cellular components, enrichment was observed in external side of plasma membrane, plasma membrane signaling receptor complex, and T cell receptor complex. Moreover, Molecular function analysis revealed enrichment in antigen binding, immune receptor activity, and cytokine receptor binding (Fig. 4D–E). GSEA analysis identified the key pathways in the high-risk group as KEGG_CALCIUM_SIGNALING_PATHWAY, KEGG_DILATED_CARDIOMYOPATHY, KEGG_ECM_RECEPTOR_INTERACTION, KEGG_FOCAL_ADHESION, and KEGG_NEUROACTIVE_LIGAND_RECEPTOR_INTERACTION. Conversely, the primary pathways in the low-risk group were identified as KEGG_GLYCEROLIPID_METABOLISM, KEGG_METABOLISM_OF_XENOBIOTICS_BY_CYTOCHROME_P45, KEGG_OXIDATIVE_PHOSPHORYLATION, KEGG_PARKINSONS_DISEASE, and KEGG_RIBOSOME. These results demonstrate that alterations in the immune microenvironment and immune functions are the primary features of the high-risk group in melanoma.

Fig. 4

Gene enrichment analysis of melanoma risk score groups. A Heatmap displays the expression levels of the top 20 differentially expressed genes in the high-risk group and low-risk group. B Volcano plot illustrates the fold change and distribution of differentially expressed genes between the high-risk group and low-risk group. C KEGG enrichment analysis demonstrates the enrichment of differentially expressed genes in KEGG pathways. D-E GO enrichment analysis demonstrates the enrichment of differentially expressed genes in GO pathways. F-G GSEA enrichment analysis demonstrates the enrichment of differentially expressed genes in GSEA pathways in different risk score groups

Melanoma is characterized by immune suppression and evasion, leading us to analyze the correlation between immune cells in melanoma and patients’ prognosis. Our findings indicate that T cell CD4 + memory activated, activated NK cells, and NK cell activated can significantly improve the prognosis of melanoma patients (Fig. S1). Conversely, T cell CD4 + naïve, activated mast cells, and myeloid dendritic cells significantly worsened the prognosis of melanoma patients (Fig. S1). These results highlight the critical correlation between immune cells in melanoma and patient prognosis, suggesting strong immune regulation as a key factor in melanoma as presented by previous studies.

Immune microenvironment analysis in different risk score groups

We further conducted a detailed analysis of the immune microenvironment in different risk score groups of melanoma. Using algorithms such as XCELL, TIMER, QUANTISEQ, MCPCOUNTER, EPIC, CIBERSORT-ABS, and CIBERSORT, we evaluated the immune cell infiltration in each melanoma sample. We found that the majority of immune cell infiltrates were reduced in samples from the high-risk group (Fig. 5A). Analysis of the immune microenvironment components demonstrated significantly decreased ESTIMATE score, ImmuneScore, and StromalScore in the high-risk group samples (Fig. 5B–D). Furthermore, almost all immune cells and immune functions were significantly suppressed in the high-risk group (Fig. 5E, F). Additionally, all immune checkpoint genes were notably downregulated in the high-risk group (Fig. 5G). These results indicate that the high-risk group of melanoma is in a state of pronounced immune suppression and exhaustion.

Fig. 5

Immune microenvironmental status in different melanoma risk score groups. A Evaluating the immune cell infiltration in each melanoma sample using algorithms such as XCELL, TIMER, QUANTISEQ, MCP COUNTER, EPIC, CIBERSORT-ABS, and CIBERSORT. B Analyzing the ESTIMATE Score in different risk score groups. C Analyzing the Immune Score in different risk score groups. D Analyzing the Stromal Score in different risk score groups. E Analyzing the immune cell content in different risk score groups. F Analyzing the degree of immune function activation in different risk score groups. G Analyzing the expression of immune checkpoint genes in different risk score groups. H Analyzing the targeted drug sensitivity in the low-risk score group. I Analyzing the targeted drug sensitivity in the high-risk score group

Drug sensitivity analysis of different risk score groups

Given the close relationship between melanoma and the immune system, targeted immunotherapy drugs and chemotherapy drugs may play a significant role. However, due to the distinct immune microenvironment in different risk score groups, it is crucial to identify appropriate treatment options for each patient. We performed drug sensitivity analysis on each melanoma sample using R. limma, R. oncoPredict, and R. parallel, with the half maximal inhibitory concentration (IC50) as the evaluation criterion for drug sensitivity. The results revealed that patients in the low-risk group exhibited greater sensitivity to Temozolomide, Ribociclib, 5-Fluorouracil, and Dasatinib (Fig. 5H), while patients in the high-risk group demonstrated better sensitivity to Selumetinib, Trametinib, Ulixertinib, and Lapatinib (Fig. 5I).

Analysis of tumor mutation burden in different risk score groups

TMB refers to the quantity of genetic mutations occurring within cancerous cells. The more gene mutations or aberrant sites present, the higher the likelihood of detection by the body’s immune system, consequently triggering the body’s immune defense mechanism. This could lead to immune cells searching for and killing cancer cells. Therefore, it can be inferred that higher TMB values may potentially lead to more favorable outcomes with the use of immune-based medications. Based on the R. ggpubr and R. reshape2 packages, we calculated the TMB for each melanoma sample and conducted a specific analysis of the mutated genes. Our results indicated that both the high-risk group and the low-risk group exhibited TMB levels at a relatively high percentage (> 90%), suggesting that overall, melanoma could benefit from immune treatment (Fig. 6A). In comparison to the high-risk group, the low-risk group demonstrated higher TMB, particularly in the TNN and MUC16 genes, suggesting that the low-risk group may be more likely to benefit from immune drugs targeting these specific gene sites (Fig. 6A–C). Furthermore, Pearson correlation analysis revealed a decrease in TMB value with increasing risk score (Fig. 6D). Additionally, patients with low TMB had significantly poorer prognosis compared to those with high TMB, possibly due to the inferior response to immune treatment in low TMB patients (Fig. 6F). Given the correlation between risk score and TMB value (Fig. 6D), we attempted to combine risk score and TMB to more accurately predict the prognosis of melanoma patients. The results showed that patients with low TMB combined with high risk scores had the worst prognosis, while patients with high TMB combined with low risk scores had the best prognosis, and those with high TMB combined with high risk scores and low TMB combined with low risk scores appeared in between the two (Fig. 6F). These results indicate an association between high risk scores and low TMB values in melanoma patients and suggest that risk score and TMB values can be used together to predict the prognosis of melanoma patients.

Fig. 6

Analysis of tumor mutation burden in different melanoma risk score groups. A Waterfall plot displays the tumor mutation burden in the high-risk score group. B Waterfall plot displays the tumor mutation burden in the low-risk score group. C Box plot displays the distribution of tumor mutation burden among different risk score groups. D Pearson analysis demonstrates the correlation between risk score and tumor mutation burden. E Survival curve displays the prognosis of melanoma patients in different tumor mutation burden groups. F Survival curve displays the prognosis of melanoma patients predicted by the combined effect of tumor mutation burden and risk score

New molecular classification of melanoma based on 2-DRL signatures

The immune microenvironment of melanoma can vary significantly, and the 2-DRL model has been proven to distinguish between high-risk and low-risk subtypes of melanoma. Therefore, we attempted to use the 2-DRL model to construct a novel molecular classification for guiding melanoma treatment. Using the R.limma package and the R. ConsensusClusterPlus package, we applied unsupervised clustering methods based on the 2-DRL model to classify melanoma samples. The cumulative distribution function (CDF) curve and triangle area plot demonstrated that the CDF curve was smoother and reached its maximum value when K was set to 3, indicating a more reliable clustering result (Fig. 7A, B). The samples of different subtypes of melanoma were also evenly distributed when K was 3 (Fig. 7C). The clustering consistency among different subtypes was good when K was 3 (Fig. 7D). Based on the Kaplan–Meier method, survival curves were plotted, revealing significant differences in patient prognosis among the three molecular subtypes of melanoma. The C2 group had the best prognosis, the C3 group had the worst prognosis, and the C1 group fell in between the two (Fig. 7E). The Sankey diagram showed that the majority of patients in the C3 group belonged to the high-risk group, while the majority of patients in the C2 group belonged to the low-risk group. The number of high-risk and low-risk patients in the C1 group was relatively comparable (Fig. 7F). PCA (Principal Component Analysis) using the R.ggplot2 package demonstrated clear differentiation between high-risk and low-risk patients, as well as among the three molecular subtypes (Fig. 7G–H).

Fig. 7

A new molecular classification of melanoma based on 2-DRL signatures. A-B CDF (cumulative distribution function) curve and triangle area plot display the reliability of unsupervised clustering results. C-D Heatmap illustrates the distribution of patients in different molecular subtypes. E Survival curve demonstrates the prognosis of patients in different molecular subtypes. F Sankey diagram shows the correspondence between patients in different molecular subtypes and different risk score groups. G PCA (Principal Component Analysis) displays the discriminatory ability of the 2-DRL model for melanoma patients with different risk scores. H PCA demonstrates the discriminatory ability of the 2-DRL model for melanoma patients with different molecular subtypes

Immune signatures of various molecular subtypes of melanoma

We further analyzed the immune characteristics of each molecular subtype of melanoma. Analysis of immune microenvironment components showed that the ESTIMATEScore, ImmuneScore, and StromalScore in the C2 group were the highest among the three groups, while the C1 group had the lowest scores, with the C3 group falling in between (Fig. 8A–C). Analysis of immune cell infiltration revealed a significant increase in immune cell infiltration in the C2 group, while both the C1 and C3 groups exhibited lower levels of immune cell infiltration (Fig. 8D). Additionally, nearly all immune checkpoint genes were upregulated in the C2 group and downregulated in the C1 and C3 groups (Fig. 5E). These findings suggest that the C2 group of melanoma is in an immune-activated state, which could explain its better prognosis. On the other hand, the C1 and C3 groups appear to be in an immune-suppressed and immune-exhausted state, potentially contributing to their poorer prognosis.

Fig. 8

Immune signatures of various molecular subtypes of melanoma. A Analyzing the ESTIMATE score in different molecular subtypes. B Analyzing the Immune Score in different molecular subtypes. C Analyzing the Stromal Score in different molecular subtypes. D Evaluating the immune cell infiltration in each melanoma sample using algorithms such as XCELL, TIMER, QUANTISEQ, MCP COUNTER, EPIC, CIBERSORT-ABS, and CIBERSORT. E Analyzing the expression of immune checkpoint genes in different molecular subtypes. F Analyzing the targeted drug sensitivity in cluster1 molecular subtype. G Analyzing the targeted drug sensitivity in cluster2 molecular subtype. H Analyzing the targeted drug sensitivity in cluster3 molecular subtype

Drug sensitivity analysis of various molecular subtypes of melanoma

We further analyzed the drug sensitivity of each molecular subtype of melanoma to guide clinically targeted treatment. The results showed that the C1 group exhibited good sensitivity to Zoledronate, UMI-77, Nilotinib, and Cytarabine (Fig. 8F). The C2 group showed good sensitivity to Ribociclib, XAV939, Topotecan, and Ruxolitinib (Fig. 8G). The C3 group demonstrated good sensitivity to VX-11e, Ulixertinib, Trametinib, and Afatinib (Fig. 8H). These findings can provide valuable references for clinical decision-making by healthcare professionals regarding treatment options.

Discussion

Melanoma, due to its strong immunogenicity, serves as an excellent model for evaluating innovative immunotherapies such as checkpoint inhibitors, anti-cancer vaccines, and engineered chimeric antigen receptor T-cell (CAR-T cell) therapy. Additionally, melanoma may be susceptible to the influence of novel checkpoint inhibitors targeting B and T lymphocyte attenuator (BTLA), T-cell immunoglobulin, and mucin domain 3 (TIM-3), and lymphocyte activation gene 3 (LAG-3). This area of research remains an ongoing field of investigation [15, 16]. Despite significant advancements in cancer immunotherapy, a substantial proportion of melanoma patients do not respond or experience relapse due to primary or acquired resistance. Treatment failure rates range from 40 to 65% among patients receiving anti-PD-1 therapy, with over 70 cases of treatment failure reported in melanoma patients receiving anti-CTLA-4 therapy [17, 18]. Recent studies have discovered a novel form of cell death known as Disulfidptosis, which is closely related to the immune microenvironment and may hold potential therapeutic value and targeted treatment options for melanoma [19]. Therefore, exploring the biological characteristics associated with Disulfidptosis-related genes could open up a new research field. We constructed a 2-DRL prognosis model based on DRG and observed favorable performance of the model in predicting melanoma prognosis, as demonstrated by survival curves and ROC curves. Furthermore, we discovered a significant decrease in ESTIMATEScore, ImmuneScore, and StromalScore in the high-risk group of melanoma, indicating reduced immune cell infiltration and a state of apparent immune suppression and exhaustion. Improving the immune response in these high-risk patients remains an unresolved challenge. Therefore, through drug sensitivity analysis, we found that low-risk patients demonstrated higher sensitivity to Temozolomide, Ribociclib, 5-Fluorouracil, and Dasatinib, while high-risk patients exhibited better sensitivity to Selumetinib, Trametinib, Ulixertinib, and Lapatinib. The 2-DRL prognosis model not only predicts melanoma patients’ prognosis but also enables the analysis of the immune microenvironment in different risk score groups to select appropriate targeted drug treatments, potentially significantly improving the cure rate of melanoma immunotherapy and prolonging patient survival.

The majority of melanoma patients are diagnosed with primary tumors, and over the past decade, the prognosis for melanoma patients has seen some improvement with the introduction of novel systemic therapies [20, 21]. However, the median survival of melanoma patients still remains around three years, which is relatively low. It is important to note that a large proportion of patients who ultimately succumb to melanoma initially have early-stage disease, indicating that a subset of patients harbors aggressive tumors [22, 23]. Therefore, the identification of optimal treatment strategies for early-stage melanoma is crucial for improving outcomes. Using univariate and multivariate COX regression analyses, we constructed a column diagram combining clinical and pathological factors with the 2-DRL prognosis model. The accuracy of the column diagram was confirmed with a C-index of 0.714. The column diagram we constructed is capable of predicting the 1-year, 3-year, and 5-year survival rates of melanoma patients, providing favorable and accurate references for clinical decision-making.

Diagnosis of cutaneous melanoma has been established over the past several decades, primarily based on histopathological features completed with a relatively simple set of immunohistochemical markers. However, this has recently undergone a profound change; widespread screening programs for cutaneous melanoma in dermatology detect precursor lesions at a higher frequency, necessitating highly sensitive molecular testing to confirm malignancy [24]. On the other hand, the genetic background of hereditary melanoma is becoming increasingly complex, once again requiring more sophisticated genomic testing [24]. Last but not least, due to cutaneous melanoma being the most metastatic human cancer to be discovered at an early stage, there is a clinical need for more precise prediction, which can be based on newly developed genetic testing. The expression features of genes have the potential to improve the prediction of the biological behavior of melanoma by objectively defining “high risk” at the molecular level [24]. Identification of genes driving cutaneous melanoma is helpful in developing targeted therapies, making predictive genetic characterization a routine practice [25, 26]. However, the most effective treatment for cutaneous melanoma is immunotherapy; unfortunately, although our understanding of the immune genomic features of cutaneous melanoma has greatly increased in the past few years, its predictive markers have not yet entered clinical practice [25, 26]. We have constructed three new molecular subtypes of melanoma through the 2-DRL prognosis model, which can more accurately reflect the immune microenvironment status of each molecular subtype, to refine the personalized precision treatment of melanoma. Immune microenvironment component analysis and immune infiltration analysis show that melanoma C2 subtype is in an immune-activated state, which may be the reason for its better prognosis, while C1 and C3 subtypes are in an immune-suppressive and immune-exhausted state, which may lead to a worse prognosis. We further analyzed the targeted drug sensitivity of each molecular subtype of melanoma to guide clinical precision treatment. The results show that the C1 subtype has better sensitivity to Zoledronate, UMI-77, Nilotinib, and Cytarabine. The C2 subtype has better sensitivity to Ribociclib, XAV939, Topotecan, and Ruxolitinib. The C3 subtype has better sensitivity to VX-11e, Ulixertinib, Trametinib, and Afatinib. These results can provide effective references for clinical treatment decisions. Although our analysis is complete and reliable, the therapeutic drugs for different subtypes obtained through analytical screening still need to be verified by further basic and clinical experiments.

Long non-coding RNAs (lncRNAs) have emerged as novel pharmacological targets for the treatment of melanoma due to their crucial roles in the disease [27, 28]. There is still a lack of comprehensive research on these molecules (USP30-AS1 and LINC02560) in melanoma. Based on the achievements of the 2-DRL prognostic model, USP30-AS1 and LINC02560 hold promise as novel therapeutic targets for further investigation in melanoma treatment. Significant progress has been made in recent years in optimizing delivery strategies and chemical modifications of ncRNA targets, including polymers, lipid nanoparticles, and extracellular vesicles [10]. These ncRNA targeting approaches are easier and more effective than designing specific protein-binding inhibitors, with antisense oligonucleotides (ASOs) and RNA interference (RNAi) being the most widely used methods[29]. Based on these research advances, we can design targeted ASOs and RNAi according to the RNA sequences of USP30-AS1 and LINC02560, and try targeted delivery through polymers, lipid nanoparticles, and extracellular vesicles.

This study still has some limitations. Nowadays, research on distinguishing molecular subtypes of various cancers is in full swing. For example, Zhang et al. developed an integrated algorithm and process to identify immune-related lncRNAs to decipher gene relationships [30]. However, due to the non-disclosure of the relevant code, we did not use the relevant new algorithm to build the prognostic model and compare with the advantages of our existing algorithm. In the future, we should try to obtain more accurate models when we can obtain or reproduce the relevant new algorithms. Another limitation of this study is the lack of external datasets and experimental validation. Since there is only one control patient in the TCGA-SKCM data, we may obtain some false positive or false negative differential genes in the differential analysis. In addition to the TCGA dataset, there is currently a lack of complete clinical data and lncRNA expression profiles related to melanoma, which makes it impossible for us to perform external validation. In the future, we should collect enough melanoma samples for lncRNA sequencing to further verify our conclusions, and we should conduct in-depth research on the genes in the model through oncology-related experiments.

In conclusion, we have established a two-DRL prognostic model and a novel molecular subtyping for melanoma, revealing the immune microenvironment status and targeted drug sensitivity of melanoma patients in different risk score groups and molecular subtypes. This provides valuable guidance for clinical treatment and identifies significant DRL targets for future in-depth research.