Disulfidptosis-Related Classification Patterns and Tumor Microenvironment Characterization in Skin Cutaneous Melanoma

Abstract

Aim: To identify distinct disulfidptosis-molecular subtypes and develop a novel prognostic signature. Methods/materials: We integrated into this study multiple SKCM transcriptomic datasets from the Cancer Genome Atlas database and Gene Expression Omnibus dataset. The consensus clustering algorithm was applied to categorize SKCM patients into different DRG subtypes. Results: Three distinct DRG subtypes were identified, which were correlated to different clinical outcomes and signaling pathways. Then, a disulfidptosis-relaed signature and nomogram were constructed, which could accurately predict the individual OS of patients with SKCM. The high-risk group was less sensitive to immunotherapy than the low-risk group. Conclusion: The signature can assist healthcare professionals in making more accurate and individualized treatment choices for patients with SKCM.

Graphical abstract

Keywords::

• disulfidptosis
• immunotherapy
• risk model
• skin cutaneous melanoma
• tumor microenvironment

Skin cutaneous melanoma (SKCM) is a malignant tumor originating from melanocytes and is the third most common skin malignancy [1]. The total number of new cases of melanoma worldwide in 2020 is estimated to be 325,000 and 57,000 deaths [2]. SKCM has become one of the major problems that endanger human health. SKCM is a complex disorder characterized by a high mutational load, extensive genetic heterogeneity, and complex tumor microenvironment (TME) interactions that place it among the most aggressive types of cancer [3,4]. The highly invasive and metastatic nature of SKCM is responsible for the 5-year overall survival (OS) of only 23% of CM patients [2]. Despite having similar tumor grading and identical pathological staging, the survival outcomes of SKCM patients can vary significantly due to distinct genetic characteristics [5]. Hence, in order to enhance the result and forecast of melanoma more effectively, it is imperative to delve deeper into the exploration of more precise and specific targets as well as molecular profiling.

SKCM is regarded as a highly aggressive form of cancer. Based on the 2020 worldwide cancer data, SKCM stands at the 19th position in terms of prevalence among the most frequent types of cancer [6]. The incidence of new cases has surged to 324,635, resulting in 57,043 fatalities [6]. In early-stage melanoma, when surgery was performed promptly, a positive outcome was observed, with survival rates reaching 95% after 10 years, compared with less than 20% for metastatic melanoma. Immunotherapy is the primary approach for managing advanced melanoma [7]. In the past few years, the field of CM therapy has experienced remarkable advancements and entered a fresh era with the development of CM immunotherapy. Immunotherapy, in contrast to traditional chemotherapy, can induce an unparalleled and enduring reaction in individuals suffering from advanced cancer. Nevertheless, this response occurs solely in a comparatively limited group of individuals, and the impact significantly differs among patients with SKCM. Researchers are motivated by these clinical obstacles to discover novel methods for predicting which patients possess inherent resistance to targeted therapy and immunotherapy. By doing so, it can provide improved guidance for the clinical management of patients and encourage the rational utilization of clinical resources.

Recently, a group of scientists identify a novel programmed cell death (PCD) mode induced by aberrant accumulation of intracellular disulfides and define it as disulfidptosis [8]. They identified several genes remarkably associated with disulfidptosis via whole-genome CRISPR-Cas9 screen, including SLC7A11 and its chaperone SLC3A2, and various components of mitochondrial oxidative phosphorylation system. Under glucose starvation, high expression of SLC7A11 in kidney cancer cells accelerates nicotinamide adenine dinucleotide phosphate (NADPH) depletion in the cytoplasm, intracellular disulfide accumulation, and ultimate disulfidptosis. Of note, the pharmacological blockade of glucose uptake by GLUT inhibitors was proved to exert cancer killing effects by the promoting disulfidptosis in SLC7A11-high tumor cells, highlighting the therapeutic utility of disulfidptosis-induction strategies in cancer treatment [8]. Disulfidptosis is correlated with prognosis, the tumor microenvironment (TME), immune evasion, and therapeutic outcomes in other forms of cancer, including colon cancer [9,10], bladder cancer [11], and renal clear cell carcinoma [12]. Li et al. [10] indicated that lack of effective immune infiltration were observed in the high-risk score group, indicating a immune-exclusion TME phenotype, with poorer survival. Low-risk score group was also linked to increased anti-tumor immune cell infiltrations and elevated activation status of anti-tumor immunity. Despite its recent emergence, the association between disulfidptosis and prognosis and TME in SKCM is unclear.

The genomic characteristics of disulfidptosis-related genes (DRGs) specific to SKCM were thoroughly investigated in this study. Unsupervised consensus clustering identified three distinct disulfidptosis expression patterns based on DRGs. We clarified the variations in prognosis, clinical traits, and immune characteristics among the three clusters. Furthermore, we investigated the predictive function of these DRGs among the three clusters, conducted functional analysis on genes that were differentially expressed among different clusters, and developed a prognostic model. This signature helps measure the disulfidptosis-related features, where a high-risk score indicates a poor prognosis and a lower tumor mutation burden (TMB) in patients with SKCM. Afterward, we examined scores for evaluating the tumor microenvironment (TME), associations between TMB, and disparities in chemotherapy sensitivity within the two risk groups. The findings indicate that DRGs have a significant impact on SKCM, aiding in the assessment of patient prognosis and their reaction to chemotherapy and immunotherapy. Additionally, these genes have the potential to serve as collaborative targets for enhancing the effectiveness of SKCM treatment.

Materials & methods

Individuals & collections of data

We examined melanoma datasets (TCGA-SKCM and GSE65904) in a pair of databases, the Cancer Genome Atlas (http://portal.gdc.cancer.gov/) and GEO (https://www.ncbi.nlm.nih.gov/geo/). The normal human skin transcriptome data stored in the Genotype-Tissue Expression (GTEx) database were downloaded from the UCSC Xena database (http://xena.ucsc.edu/) and used as a control (TPM format). Samples without substantial clinicopathological or survival data were eliminated from subsequent analysis. Finally, the TCGA-SKCM dataset (468 samples) and GSE65094 dataset (210 samples) were selected as the training and validation sets, respectively.

Unsupervised clustering analysis of DRGs

After excluding unexpressed DRGs in certain samples, we selected a total of 9 DRGs for model construction, which were obtained from previous studies [8]. Consensus Clustering, a widely used research technique for classifying cancer subtypes, is an unsupervised clustering approach. By utilizing various omics data sets, it can classify samples into distinct subcategories, enabling the identification of novel disease subtypes or facilitating comparative analysis among different subtypes. To differentiate various molecular subtypes, the “ConsensusClusterPlus” R package was employed for performing consensus clustering according to the expression level of 9 DRGs. Furthermore, to further examine the distinctions among subcategories. To investigate the distribution of various subtypes, we utilized the t-distributed stochastic neighbor embedding (tSNE) technique and employed the R package “Rtsne” and “umap” to assess the impact of classification. Moreover, we examined the level of immune cell infiltration among various subcategories. Different subtypes were analyzed using the R package called “heatmap” to examine the expression levels of DRGs, age, gender, tumor location, stage, Breslow depth, ulceration, and tumor status. In the end, the analysis of the enriched pathways between various subtypes was performed using the “GSVA” R package and visualized as heatmaps.

Establishment & validation of DRG-related signature

To identify the different subtypes’ differentially expressed genes, the R package called “limma” was utilized using the criteria of FDR <0.05 and logFC >0.585. Following the acquisition of the distinctively expressed genes among each subtype, we selected the common genes for further analysis. The samples in TCGA database were used as the training set, and the samples in GSE65904 dataset were used as the testing set. To identify genes associated with prognosis (p < 0.05), the Univariate Cox regression analyses were utilized to discover intersecting genes in the training cohort. The “glmnet” R package was utilized to conduct Least Absolute Shrinkage and Selection Operator (LASSO) Cox regression analysis using differentially expressed genes associated with prognosis in the training set. The risk score was calculated by multiplying the sum of the expression levels for each gene with the LASSO regression coefficient for each gene. The samples were categorized into high and low groups based on the median DRG-score. Afterward, we examined the AUC of the training set and the validation set.

Clinical significance of the DRG-related signature

Univariate and multivariate Cox regression analysis was conducted to validate the independence of the DRG-score as a prognostic predictor, considering both the risk score and clinicopathological variables. The findings were unveiled on the map of the woodland. Next, we performed a classification analysis to investigate if the DRG-score maintains its predictive accuracy in different subgroups defined by various clinical variables.

Development & verification of a nomogram scoring model

Based on the autonomous prognosis result, a prognostic nomogram was generated using the clinical features and the DRG-score with the assistance of the “rms” package in R. Each variable in the nomogram is assigned a specific score, and the overall score is calculated by summing up the scores of all variables for each sample. ROC curves were used to evaluate the nomogram’s performance in predicting survival rates at 1, 3, and 5 years. The calibration plots of the nomogram were utilized to depict the prognostic significance of the expected survival events at 1, 3, and 5 years in comparison to the observed outcomes.

Analysis of tumor microenvironment

Immune cell infiltration analysis was conducted using the ssGSEA algorithms. Wilcoxon signed-rank test was used to analyze the distinct composition of immune infiltrating cells in the high- and low-risk groups. We examined the relationship between immune cells and 8 crucial genes by conducting a correlation analysis. Simultaneously, we examined the association between the two predictive risk categories and the TME. A box plot was generated using the “estimate” package to compare the TME scores between the different DRG-score groups.

Anticipating the outcome of immunotherapy

The Tumor Immune Dysfunction and Exclusion (TIDE) method was utilized to determine the anticipated reactions to the immune checkpoint inhibitors by comparing the variations in TIDE scores for each sample between the distinct groups. In each sample from the TCIA database, an immunophenoscore (IPS) was created, which is a notable indicator of response to anti-PD-1 and anti-CTLA-4. Subsequently, the IPS was compared across risk groups in TCGA-SKCM to investigate the correlation between the risk score and IPS. A box plot was generated using the “reshape2” and “ggplot2” packages to compare the expression levels of immune checkpoints between the high-risk and low-risk groups. Additionally, the correlation between immune cells and the signature genes was analyzed.

Evaluation of drug responsiveness

Patients were divided into two subgroups based on their DRG-score, and the prediction of drug sensitivity was made for different medications. Drug prediction was performed using the “pRRophetic” R package. The Wilcoxon signed-rank test was employed to investigate the variation in IC50 values among various risk categories. The analysis was conducted using the R package “ggplot2”.

The relationship between the DRG-score & TMB

Data on somatic mutation data of SKCM was obtained from the TCGA database for the analysis of gene mutation. The gene mutation in various risk subgroups was analyzed using the “Maftools” R package. Subsequently, the analysis of the relationship between the DRG-score and the tumor mutation burden (TMB) was performed. Furthermore, we performed survival analysis among various TMB subcategories to investigate the influence of TMB status on the prognosis of patients with SKCM. Afterward, we merged TMB and DRG-score to conduct survival analysis on SKCM patients.

Detection of differentially expressed genes & analysis of functional enrichment

To detect differentially expressed genes (DEGs) in the different risk subgroups, we employed the “limma” software package, applying the following criteria: |Fold Change| >1.5 and false discovery rate (FDR) <0.05. The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis was conducted using the “clusterPrfiler” with a significance threshold of FDR <0.05.

Statistical analysis

The data analysis was conducted using R software and R Bioconductor packages. The version is 4.3.0. The distinctions between non-parametric and parametric methods were evaluated using the Wilcoxon test, Kruskal-Wallis test, t-test, or one-way ANOVA. The ROC curve was used to confirm the accuracy of the model. The prognosis of the survival curve was analyzed using the Kaplan-Meier test and Log-rank test, which were employed to evaluate distinctions among the groups. Statistical analyses were conducted on both sides, and a significance level of p < 0.05 was used.

Results

Development of DRG molecular subtypes of melanoma

We employed the unsupervised clustering method to detect distinct regulatory patterns of SKCM by analyzing the expression levels of DRGs. After analyzing 468 samples, it was determined that the ideal number of clusters is 3 in the cluster analysis. At K = 3, the groups exhibited the least disparity, while the dissimilarity beyond the group was the most significant. As a result, we successfully classified melanoma patients into three distinct subcategories, specifically A, B, and C (Figure 1A). The stable distribution of melanoma patients was closely indicated by the relative change in the area under the CDF curve when dividing them into 3 subtypes (Figure 1B, C). The tSNE, UMAP, and PCA analysis indicated that the A, B, and C subtypes can be differentiated from one another. Figure 1D–F demonstrates that our classification using DRGs for subtype analysis exhibits superior typing capability.

Figure 1. Using unsupervised clustering analysis to cluster SKCM based on differentially expressed genes.

(A) The heatmap of the consensus matrix delineates three clusters (k = 3) and their corresponding correlation region. (B) a CDF curve of K = 2–10. (C) The relative change in area under the CDF curve of K = 2–10. (D–F) PCA, tSNE, and UMAP identified three DRG clusters.

CDF: Cumulative distribution function; PCA: Principal component analysis; tSNE: t-distributed stochastic neighbor embedding; UMAP: Uniform manifold approximation and projection.

Comparative examination of the three DRG subtypes

The Kaplan-Meier analysis showed that the C subtype had the highest prognosis, followed by the B subtype, while the A subtype had the worst survival outcome (Figure 2A). In a heat map, we display the levels of DRGs and clinical characteristics, including age, gender, tumor location, stage, Breslow depth, ulceration, and tumor status for the A, B, and C subtypes. The A subtype exhibited the highest level of DRG expression, followed by the B subtype, and the C subtype showed the lowest level (Figure 2B). Next, the disparities in biological behavior among different subgypes were examined through gene set variation analysis (GSVA) (Figure 2C–E). Upon comparing the A and B subtypes, it was observed that the B subtype exhibited a considerably higher enrichment in apoptosis, NOD-like receptor signaling pathway, and Protein export than the A subtype. The A subtype exhibited considerably greater enrichment in the Aminoacyl-tRNA biosynthesis, base excision repair, and RNA polymerase compared with the C subtype when comparing the A and C subtypes. Pathways like the T cell receptor signaling pathway, Natural killer cell-mediated cytotoxicity, and chemokine signaling pathway showed a slightly higher enrichment of the C subtype compared with the A subtype. In pathways Intestinal immune network for IgA production, Antigen processing and presentation, and Primary immunodeficiency, the C subtype exhibits a higher enrichment level compared with the B subtype when comparing the B subtype and C subtype. The B subtype showed a higher level of enrichment in pathways related to the chemokine signaling pathway, pentose phosphate pathway, and the citric acid cycle compared with the C subtype.

Figure 2. Differentially expressed gene-based molecular clusters with distinct clinical features and TIME landscapes.

(A) Survival analysis was conducted on three DRG clusters. (B) Expression profiles of DRGs and clinicopathological characteristics between clusters. (C–E) The heat map of GSVA revealed the variations in pathways among the three clusters. (F) Abundances of infiltrating immune cells between three clusters.

Furthermore, we examined the three SKCM subtypes to investigate the 23 infiltrating immune cell categories (Figure 2F). The outcome indicated that the majority of the invading immune cells exhibited notable variations between the three DRG subtypes, excluding CD56 dim NK cells and monocytes. In terms of immune cell infiltration, the C subtype exhibited the highest level while the A subtype had the lowest level (Figure 2F).

Establishment & validation of DRG-related signature

DEGs associated with DRG subtypes were utilized to create the DRG-related signature. We examined the genes that were expressed differently among the subtypes A, B, and C. Between the A and B subtypes, there were 504 genes expressed differentially. Additionally, there were 1626 genes expressed differentially between the A and C subtypes. Furthermore, there exists a total of 463 genes expressed differentially between the B and C subtypes. By intersecting the sets of differentially expressed genes from the three subtypes, we identified a total of 37 genes that were co-expressed across all three subtypes (Figure 3A). Next, we conducted a uniCox analysis to assess the survival importance of these genes, resulting in the identification of 26 genes meeting the criteria of p < 0.05 (Figure 3B). Afterward, we conducted LASSO analysis on 26 DRG subtype-associated prognostic DEGs to establish a predictive model (DRG-score), resulting in the identification of eight noteworthy genes (MGAT5B, SEMA6A, OCA2, SRGN, GPR143, FGL2, PTPRC, and CDH3) (Figure 3C & D). These genes were utilized from the training set data to calculate the DRG-score, with the following formula: Risk score = expMGAT5B × 0.00068 + expSEMA6A × 0.05199 + expOCA2 × 0.03466 + expSRGN × (-0.12988) + expGPR143 × 0.03095 + expFGL2 × (-0.01660) + expPTPRC × (-0.03695) + expCDH3 × 0.01336. The training set consisted of 454 patients who were stratified into high and low DRG-score groups based on the score formula’s median score calculation. The distribution of patients in the three DRG clusters and two DRG-score groups was shown in Figure 3E. To ascertain the general distribution of melanoma samples in subpopulations with low and high risk, PCA and t-SNE analyses were performed. It is possible to effectively distinguish the patients in the two subgroups (Figure 3F & G). The analysis of survival demonstrated that the samples in the two subgroups have distinct survival statuses. The OS rates were found to be lower in the high-risk subgroup (p < 0.05; Figure 3H). For 1-, 3-, and 5-year survival, the training set’s time-dependent ROC areas are 0.723, 0.728, and 0.701 (Figure 3I).

Figure 3. Development of the differentially expressed gene-related signature in the training set.

(A) The Venn diagram displays the overlap of differentially expressed DRGs. (B) The Univariate Cox regression analysis of 26 genes is represented in the forest plot. (C, D) The prognostic genes were analyzed using LASSO-Cox regression and the partial likelihood deviance. (E) The distribution of patients in the three DRG clusters and two DRG-score groups. (F, G) Principal Component Analysis (PCA) and t-distributed stochastic neighbor embedding (tSNE) reveal a clear distinction in transcriptomes among the two subcategories. (H) Analysis of overall survival (OS) in TCGA cohort using Kaplan-Meier method. (I) Receiver Operating Characteristic curves for predicting overall survival at 1 year, 3 years, and 5 years.

Afterward, the predictive model was utilized on the validation set. The PCA and t-SNE analyses demonstrated a distinct separation between the two groups at risk (Figure 4A & B). Figure 4C demonstrates that in the validation set, the high-risk group exhibited a significantly diminished prognostic value compared with the low-risk group. Figure 4D displays the time-dependent ROC areas for 1-, 3-, and 5-year survival in the validation set, which are 0.719, 0.753, and 0.717, respectively. To further confirm the accuracy of our signature, we compared the C-index, RMS, and AUC values of our signature with those of eight published risk models. Compared with the C-index and RMS of eight published risk models, our signature has great advantages (Supplementary Figure 1A & B). Furthermore, the AUC values of the signature were assessed in comparison to those of previously published signatures. Our signature’s AUC values for the ROC curve were superior to the published signatures in terms of the highest values (Supplementary Figure 1C–J).

Figure 4. Validation of the differentially expressed gene-related signature in the GSE65904 set.

(A, B) Principal component analysis (PCA) and t-distributed stochastic neighbor embedding (tSNE) plot of the GSE65904 set. (C) Analysis of overall survival (OS) in GSE65904 set using the Kaplan-Meier method. (D) Receiver Operating Characteristic curves for predicting overall survival at 1 year, 3 years, and 5 years.

Clinical significance of the DRG-related signature

The prognosis of SKCM was significantly associated with age, DRG-score, and stage, as shown by the analysis of univariate Cox regression (Figure 5A). After additional multivariate Cox regression analysis, it is demonstrated that the DRG-score remains an independent prognostic indicator even after accounting for other clinical traits (Figure 5B). Furthermore, to investigate the predictive importance of DRG-scores in patients with SKCM, the patients were categorized into various subgroups according to clinical factors. In general, the survival of high-risk patients was typically worse when compared with low-risk patients (Supplementary Figure 2).

Figure 5. Developing a nomogram to forecast the overall survival (OS) of individuals diagnosed with SKCM in the TCGA dataset.

(A & B) Univariate and multivariate Cox regression analysis was conducted to validate the independence of the DRG-score as a prognostic predictor. (C) A nomogram utilizing the DRG-score, along with age, tumor stage. (D) Time-based receiver operating characteristic (ROC) curves of the nomogram for the prediction of 1-, 3-, and 5-year OS. (E) Calibration plots are used for internally validating the nomogram.

Creating a nomogram for forecasting survival

To predict the 1-, 3-, and 5-year OS rates in patients with SKCM, a nomogram was created that combines the DRG-score with clinicopathological characteristics, considering the limited clinical usefulness of the DRG-score in predicting OS (Figure 5C). Our analysis of the nomogram model showed high accuracy for OS at 1, 3, and 5 years (Figure 5D). The proposed nomogram showed a similar performance to an ideal model based on the calibration plots (Figure 5E).

Evaluation of TME in different risk subpopulations

Based on the ssGSEA algorithm, the TCGA database’s gene expression matrix for SKCM was used to calculate the enrichment scores of 16 immune cells and the activity of 13 immune-related pathways. Afterward, we investigated the makeup of immune cells and immune-related pathways in various risk subcategories. The findings indicated that individuals in the high-risk category exhibited notably low levels of immune cell infiltration (Figure 6A). Similarly, all the immune-related pathways were less active in the high-risk group than in the low-risk group (Figure 6B). Furthermore, there were significant variations in the tumor immune microenvironment between low- and high-risk subpopulations. Patients classified as high-risk displayed reduced ImmuneScore, StromalScore, and ESTIMATEScore levels but higher TumorPurity levels in comparison to low-risk patients (Figure 6C–F). Additionally, we examine the correlation between immune cell enrichments and the eight genes in the prognostic signature. Based on Figure 6G, we reached the determination that the majority of immune cells exhibited a strong correlation with the chosen genes.

Figure 6. Evaluation of TME in different risk subpopulations.

(A & B) Comparison of the enrichment levels of immune-related cells and pathways in the low- and high-risk groups. (C–F) A box plot illustrates the variances in the TME score when comparing the two DRG-score groups. (G) The relationship between the model genes and immune-related cells.

Assessment of the immunotherapeutic reaction

To assess the response to immunotherapy, we conducted an analysis called tumor immune dysfunction and exclusion (TIDE). By utilizing two primary tumor immune evasion mechanisms, TIDE can forecast the response to immunotherapy. These mechanisms include the inhibition of T cell dysfunction and T cell infiltration in tumors that have low CTL levels. A greater TIDE score indicated a greater likelihood of immune evasion, suggesting that patients were less likely to experience benefits from ICI treatment. Our findings indicate that patients in the non-responder group had a risk score that was comparatively higher than the responder group, suggesting that ICI therapy may be more advantageous for DRG-low patients (Figure 7A). Utilizing the TCIA repository, we employed it to produce IPS for every SKCM specimen. Notably, patients classified as low-risk exhibited higher IPS for anti-CTLA-4, anti-PD-1, and anti-(CTLA-4 plus PD-1) compared with high-risk patients, implying more favorable immunotherapy results among individuals possessing a lower risk score (Figure 7B–D). Furthermore, the group with low risk exhibited increased expression of significant immune checkpoint biomarkers (PD-1, PD-L1, and CTLA-4) in contrast to the high-risk group (Figure 7E). The results suggest that individuals with a low-risk profile may experience greater advantages from immunotherapy.

Figure 7. Assessment of the responsiveness to immunotherapy using differentially expressed gene score.

(A) Comparison of the tumor immune dysfunction and exclusion (TIDE) in the two DRG-score groups. (B–D) Comparison of the immunophenoscore (IPS) in the two DRG-score groups, categorized by anti-PD-1, anti-CTLA-4, or anti-(CTLA-4 plus PD-1). (E) Comparison of immune checkpoint gene expression in the two DRG-score groups.

Identification of potential medications for the treatment of melanoma

In the training cohort, we attempted to establish connections between various risk categories and the efficacy of chemotherapy in the treatment of SKCM, excluding the evaluation of ICI therapy. We demonstrated that a lower half inhibitory concentration (IC50) of chemo-therapeutics like 5-Fluorouracil, Cisplatin, Dasatinib, Gemcitabine, and Ribociclib (p < 0.05) was linked to low risk, while a low IC50 such as Tamoxifen, Dihydrorotenone, and Lapatinib (p < 0.05) was associated with high risk. Hence, it is demonstrated that the DRG-score functioned as a possible indicator for chemo-responsiveness (Supplementary Figure 3).

The relationship between the DRG-score & TMB

To gain a deeper understanding of the immunological characteristics in various risk subcategories, we examined the genetic mutations. In the high (Figure 8A) and low DRG-score groups (Figure 8B), we have identified the 20 genes exhibiting the most elevated mutation rates. The findings demonstrated that the occurrence of missense mutation and Multi Hit was the predominant type of mutation. Not only were the mutation frequencies of TTN, MUC16, DNAH5, BRAF, and PCLO higher than 40% in both groups, but they were also the most prevalent mutations in both groups. Furthermore, we examined the correlation between the risk score and TMB. The level of TMB exhibited a notable increase in the subgroup with low risk compared with the subgroup with high risk (Figure 8C). To investigate the influence of TMB status on the prognosis of patients with SKCM, we additionally performed survival analysis among various TMB subcategories. Patients with a high TMB exhibited a more favorable prognosis compared with patients with a low TMB (Figure 8D). Afterward, we merged TMB and DRG-score to conduct survival analysis on patients with SKCM, and the DRG-score nullified the prognostic advantage in the high-TMB category (Figure 8E).

Figure 8. Correlation between the differentially expressed gene-score and tumor mutational burden.

(A & B) Waterfall diagram illustrating the somatic mutation of tumors in both the high-risk and low-risk groups. (C) The tumor mutational burdenvaries among the subgroups with different DRG score. (D) Survival analysis was conducted on the various groups, categorized based on tumor mutational burden. (E) Analyzing the survival of different groups categorized by both tumor mutational burden and DRG score.

Detection of DEGs & analysis of functional enrichment

To better understand the physiological functions and pathways linked to risk scores, GO and KEGG enrichment analysis was performed on DEGs between the two DRG-score groups. A total of 1448 differentially expressed genes (DEGs) were identified using the criteria of |Fold Change| >1.5 and FDR <0.05 (Figure 9A, B). The GO terms include Biological Process (BP) categories like immune receptor activity, Cellular Component (CC) categories like T cell receptor complex, and Molecular Function (MF) categories like production of molecular mediators of immune response (Figure 9C). Furthermore, the KEGG results revealed that the DEGs were predominantly enriched in pathways related to the immune system, such as the T cell receptor signaling pathway, natural killer cell-mediated cytotoxicity, and antigen processing and presentation (Figure 9D). It was unexpected to discover that immunity was linked to both GO terms and KEGG pathways.

Figure 9. Analysis of functional enrichment for differentially expressed genes.

(A & B) The heatmap and Volcano plot of differentially expressed genes (DEGs). (C) The analysis of GO enrichment for BP, CC, and MF terms revealed the potential role of the DEGs. (D) The possible pathways were identified through Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis.

Discussion

Melanoma, one of the most lethal forms of skin cancer, is distinguished by its significant invasiveness and mortality rates [13,14]. The treatment approaches for advanced melanoma have undergone significant transformations over the last ten years as a result of the implementation of targeted therapy and immunotherapy [15,16]. Nevertheless, even though the implementation of focused therapy and immunological checkpoint treatments exhibits considerable potential, the emergence of swift resistance continues to pose a predominantly unconquerable obstacle [17–19]. Therefore, it is imperative to investigate novel approaches for treatment, and disulfidptosis may emerge as a promising area of future research. Extensive research indicates that addressing disulfidptosis could emerge as a novel approach for treating SKCM. The regulation of the redox state is closely linked to the metabolism of disulfides, serving as a significant mechanism in the occurrence and development of tumors [20]. For instance, certain cancerous cells enhance their ability to survive and endure by modifying the redox conditions within their cells [21]. Furthermore, certain research has indicated that disulfides could hold promise in the field of cancer therapy. For instance, the anti-cancer properties of certain drugs like cisplatin and paclitaxel are manifested through their interaction with intracellular disulfides [22,23]. Moreover, extensive studies have demonstrated that overexpression of SLC7A11 in cancer suppresses ferroptosis and thus plays a crucial role in promoting tumor growth [24–26]. The discovery of disulfidptosis subverts this traditional thinking, that is, SLC7A11 also plays an essential role in promoting disulfidptosis. Because SLC7A11 not only inhibits the occurrence of ferroptosis but also promotes disulfidptosis. The treatment that promotes ferroptosis by interfering with SLC7A11 may inhibit the occurrence of disulfidptosis. Hence, maintaining a balance between ferroptosis and disulfidptosis could emerge as a novel approach for enhancing the rate of treatment response and survival in individuals with SKCM. At the same time, preclinical findings suggest that metabolic therapy using glucose transporter (GLUT) inhibitors can trigger disulfidptosis and inhibit cancer growth [8]. Future studies are required to explore the anticancer effects of disulfidoptosis inhibitors in combination with other therapies, such as immunotherapy.

SKCM exhibits significant heterogeneity, leading to diverse clinical outcomes and varying responses to treatment. To tackle this problem, we explored the possibility of disulfidptosis in impeding the advancement of SKCM and devised an innovative disulfidptosis-associated indicator to facilitate risk categorization and personalized treatment prognosis. The use of consensus clustering algorithms enables efficient analysis and identification of patient clusters with diverse characteristics in extensive datasets. Hence, we employed this unguided algorithm to detect three separate molecular subcategories (A, B, and C) by analyzing the levels of expression of 9 DRGs. Compared with subtypes A and B, subtype C showed the most favorable survival outcomes with higher levels of immune infiltration. Additionally, we employed GSVA enrichment analysis to explore the differences in biological behavior among these three subtypes. The enrichment of pathways related to immunity was observed in subtype C compared with the A and B subtypes. Afterward, we acquired a total of 37 genes that were differentially expressed and exhibited co-expression in all three subtypes. The LASSO algorithm was employed to identify eight essential genes (MGAT5B, SEMA6A, OCA2, SRGN, GPR143, FGL2, PTPRC, and CDH3) for model development and validation, after acquiring 26 genes associated with prognosis. Our risk assessment model can differentiate between groups with a high level of risk and those with a low level of risk. Our model demonstrated improved accuracy in predicting the prognosis of melanoma patients, as evidenced by Kaplan-Meier analysis, ROC curve, nomogram, and calibration plot. Afterward, we examined scores for evaluating the TME, associations between TMB, and disparities in immunotherapy and chemotherapy sensitivity within the two DRG-score groups. The findings indicate that DRGs have a significant impact on SKCM, aiding in the assessment of patient prognosis and their reaction to chemotherapy and immunotherapy. Additionally, these genes have the potential to serve as collaborative targets for enhancing the effectiveness of SKCM treatment.

Despite being commonly regarded as a cancerous malignancy with limited treatment options, the outcomes of patients have significantly improved thanks to innovative therapies that focus on vulnerable genes and immune checkpoints. This improvement is attributed to enhanced biological knowledge and groundbreaking advancements. In recent years, the employment of ICI treatment has led to a notable rise in the 5-year survival rate for melanoma patients, increasing from less than 5% to approximately 30% [27,28]. The goal of ICIs is to address the malfunctioning immune system and stimulate CD8-positive T cells to eliminate cancerous cells [29]. The current treatments have completely transformed the level of care for patients with SKCM, however, the limited rate of response and the unavoidable development of treatment resistance might hinder any future advancements in treatment results. Our research may provide additional insights into the management of melanoma. To assess the response of ICIs, TIDE, and IPS scores were calculated. The TIDE score is linked to two distinct methods of immune evasion, specifically, the impairment of cytotoxic T lymphocytes (CTLs) that infiltrate the tumor and the exclusion of CTLs. TIDE scores demonstrate the correlation between the likelihood of immune evasion by tumors and the effectiveness of ICI therapy. Based on our analysis, we noticed that SKCM patients with reduced DRG_scores exhibited decreased TIDE scores and exhibited a favorable response to anti-PD1 and anti-CTLA-4 treatment. Additionally, we examined the manifestation of various significant genes related to immunological checkpoints in both the high-risk and low-risk patient groups. The results showed that individuals with a low-risk status exhibited increased levels of immunological checkpoint gene expression, suggesting that patients with low-risk scores could potentially gain advantages from ICI treatment. Therefore, we concluded that the levels of immunological checkpoint gene expression could serve as a reliable indicator for evaluating the impact of immunotherapy in patients with SKCM. Moreover, a significant biomarker in immune checkpoint inhibitor therapy is the high tumor mutation burden (TMB-H), which helps identify tumor patients who could potentially benefit from the therapy. The basic premise is that by increasing the number of mutant proteins, it is possible to produce antigenic peptides that could potentially boost immunogenicity. The research demonstrated that patients with a low risk displayed a high TMB and a favorable prognosis. The findings showcased the precision of our predictive model in evaluating the risk of patients from an alternative standpoint.

Gene mutations are abundant during the formation and progression of tumor tissue. Altered genetic material can produce tumor antigens, which can be identified by the immune system as foreign tissues, thereby stimulating a response from immune cells [30]. The treatment of tumors greatly benefits from immunotherapy as it exploits the ability of immune cells to identify and eradicate cancer cells [31,32]. Nonetheless, tumors efficiently inhibit immune responses (immune evasion) through the activation of negative regulatory pathways related to immune balance or by acquiring characteristics that enable them to actively avoid recognition [33,34]. Highly effective immunotherapy drugs have received approval in preclinical and clinical phase I–III trials for melanoma that is aggressive, refractory, advanced, and metastatic [35]. Clinical trials are currently evaluating the efficacy of nivolumab and pembrolizumab, which are anti-PD-1 monoclonal antibodies, along with ipilimumab, an anti-CTLA-4 antibody, for the treatment of melanoma [36]. Research has indicated that frequently employed immune checkpoint inhibitors (ICIs) can enhance the survival rates of individuals with melanoma, both in terms of progression-free survival and overall survival [37,38]. Our study revealed that the high-risk group exhibited a notably reduced level of immune cell infiltration compared with the low-risk group, as indicated by the risk score model. Notably, the high-risk group had a considerably shorter survival time compared with the low-risk group. In previous investigations on melanoma, it was observed that the high immune score group exhibited a notably longer survival time compared with the low immune score group [39,40]. In this study, we put forward a hypothesis: Does the level of immune cell infiltration in melanoma have a positive correlation with patients’ survival time? To validate this issue, additional clinical data or experiments are required.

Based on the clinical trial, the literature indicated that immune therapy yielded positive results in individuals with SKCM before achieving disease control through standard chemotherapy. Our objective was to determine whether the combination of chemotherapy and immune therapy in SKCM exhibited improved efficacy, warranting further investigation. Hence, we investigated the responsiveness of different medications in individuals belonging to two distinct risk categories. The results of our research indicated that the low-risk category exhibited a considerable likelihood of responding positively to ICI treatment. Additionally, we discovered a strong correlation between low risk and the effectiveness of specific medications such as 5-Fluorouracil, Cisplatin, Dasatinib, Gemcitabine, and Ribociclib. This implies that future research can concentrate on integrated therapy for patients with SKCM. The DRG-score’s predictive model identified these medications as having the ability to potentially treat SKCM in certain circumstances.

While our prognostic signature for disulfidptosis has a remarkable capacity to detect the immune landscape and forecast prognosis in patients, it is important to acknowledge certain limitations and develop suitable approaches to address them in future studies. Analyzing data from public databases may result in biased predictions when compared with real-world scenarios. Despite our efforts to reduce this situation, additional data from patients with SKCM are required to validate the feasibility of the model and the precision of immunotherapy prognosis. Furthermore, further investigations and fundamental research are required to enhance the pertinent specifics of this study.

Conclusion

Consensus clustering was employed to detect three disulfidptosis-molecular subtypes in SKCM that exhibited distinct OS. Based on disulfidptosis-DEGs, we developed a prognostic signature that more accurately predicted patient survival outcomes. Furthermore, we tentatively established a connection between our risk model and the immune landscape. Our study findings offer valuable information for forecasting the outcomes of individuals with SKCM, and could potentially assist in their clinical treatment.

Summary points

• The discovery of disulfidptosis provides a new avenue for anticancer treatments that target the pathophysiological role of disulfide stress.

• Little is known about the impact of disulfidptosis-related genes (DRGs) on the prognosis of skin cutaneous melanoma (SKCM) and tumor microenvironment (TME).

• We integrated into this study multiple SKCM transcriptomic datasets from the Cancer Genome Atlas database and Gene Expression Omnibus dataset.

• Three distinct DRG subtypes were identified, which were correlated to different clinical outcomes and signaling pathways.

• The least absolute shrinkage and selection operator (LASSO) and Cox regression analyses identified eight genes for the construction of a prognostic signature related to disulfidptosis.

• The novel signature displayed outstanding durability and prognostic capability for the overall lifespan of individuals.

• By combining risk scores with clinicopathological features, a nomogram was constructed, which could accurately predict the individual OS of patients with SKCM.

• The high-risk group had lower recruitment of tumor growth-associated immune cells, and the anti-immune checkpoint immunotherapy was less efficacious overall than that of the low-risk group.

• SKCM patients with different risk scores have different sensitivities to antitumor drug treatment.

Author contributions

L Yang, Z-j Cao and J Tian conceived the study. L Yang, Z-j Cao and Y Zhang contributed to sample collection and patient data management. L Yang, Z-j Cao, Y Zhang and J-k Zhou analyzed the samples and data. L Yang, Z-j Cao and J Tian wrote the manuscript draft. All authors contributed to the article and approved the submitted version.

Competing interests disclosure

The authors have no competing interests or relevant affiliations with any organization or entity with the subject matter or materials discussed in the manuscript. This includes employment, consultancies, honoraria, stock ownership or options, expert testimony, grants or patents received or pending, or royalties.

Writing disclosure

No writing assistance was utilized in the production of this manuscript.

Ethical conduct of research

The authors state that they have obtained appropriate institutional review board approval or have followed the principles outlined in the Declaration of Helsinki for all human or animal experimental investigations.

Supplementary data

To view the supplementary data that accompany this paper please visit the journal website at: www.futuremedicine.com/doi/suppl/10.2217/mmt-2023-0006

Financial disclosure

This study was supported by the Natural Science Foundation of Shaanxi Province (No. 2021JQ-919), Shaanxi Provincial People’s Hospital 2020 Science and Technology Development Incubation foundation (No.2020YXM-05). The authors have no other relevant affiliations or financial involvement with any organization or entity with a financial interest in or financial conflict with the subject matter or materials discussed in the manuscript apart from those disclosed.

Additional information

Funding

This study was supported by the Natural Science Foundation of Shaanxi Province (No. 2021JQ-919), Shaanxi Provincial People’s Hospital 2020 Science and Technology Development Incubation foundation (No.2020YXM-05). The authors have no other relevant affiliations or financial involvement with any organization or entity with a financial interest in or financial conflict with the subject matter or materials discussed in the manuscript apart from those disclosed.