Figure S1 depicts a flowchart of the current work’s procedure. The Gene Expression Omnibus, The Cancer Genome Atlas (TCGA), and Genotype-Tissue Expression (GTEx) databases were used to provide gene expression (fragments per kilobase million, FPKM) and significant prognostic and clinicopathological data for SKCM. For the following analyses, four GEO melanoma cohorts (GSE19234, GSE65904, GSE78220, and GSE133713) and TCGA cohorts were collected. We retrieved the raw “CELL” files and adjusted the background and quantile standardization. The FPKM numbers of TCGA-Skin Cutaneous Melanoma (SKCM) were converted into transcripts per kilobase million (TPM) and were thought to be equivalent to those from microarrays (8).

The “Combat” technique was used to minimize batch effects after combining five datasets. We eliminated data from patients who did not have an OS; hence, 781 SKCM patients were included in the ensuing analyses. Table S1 contains detailed information on these 781 SCKM patients. Age, gender, tumor site, TNM stage, follow-up period, and survival status were all clinical factors.

Ten CRGs were found in earlier papers (7). Table S2 has complete information on these genes. To categorize patients into discrete molecular subgroups based on CRG expression, the R package “ConsensusClusterPlus” was used for consensus unsupervised clustering analysis. This clustering was done using the following criteria: First, the cumulative distribution function (CDF) curve progressively and gently grew. Second, there were no groups with small sample sizes. Finally, clustering boosted intra-group correlation while decreasing inter-group correlation. To study changes in CRGs in biological processes, gene set variation analysis (GSVA) was done using the MSigDB hallmark gene set (c2.cp.kegg.v7.2).

We investigated the correlations between genetic subgroups, clinical and pathological features, and prognosis to assess the clinical utility of the subgroups found by consensus clustering. Age, gender, tumor site, and TNM stage were among the patient characteristics. Moreover, the variations in OS across subgroups were examined using Kaplan-Meier curves built by the R programs with “survival” and “survminer” packages.

Risk Score=

where, β_i is the parameter, and gene_i is an expression of each gene.

Each patient’s immunological and stromal scores were evaluated using the ESTIMATE methodology. Furthermore, the CIBERSORT method was used to compute the percentages of 22 human immune cell types in each melanoma sample (9). A single-sample gene set enrichment analysis (ssGSEA) technique was also used to estimate the degree of immune cell infiltration in the SKCM TME (10). We also looked at the relationships between different subgroups and CLTA4 and PD-L1 expression.

DEGs between Cuprotosis subgroups were found using the R package “EdgeR” with a fold-change of 2.0 and an adjusted p-value of<0.0001. Functional enrichment studies were performed on the DEGs using the “clusterprofiler” package in R to further investigate the probable functions of cuprotosis pattern-associated DEGs and uncover related gene functions and enriched pathways.

To measure the Cuprotosis patterns of the individual tumors, the Cuprotosis score was established. The DEGs were first applied to univariate Cox regression assessment to find those associated with SKCM OS. Second, using an unsupervised clustering approach based on the expression of prognostic CRGs, the patients were divided into three subgroup groups (Cuprotosis gene subgroup A, Cuprotosis gene subgroup B, and Cuprotosis gene subgroup C) for further study. Finally, all SCKM patients were randomly assigned to training (n = 452) and testing (n = 329) sets in a 70:30 ratio, and the former was utilized to calculate the Cuprotosis-related predictive CRG score. Although, using the R package “glmnet,” we built LASSO Cox regression models based on the expression status data of cuprotosis-related prognostic genes. The least absolute shrinkage and selection operator (LASSO) is a well-known approach for assessing survival data, and it is particularly effective for studying gene expression profiles with high dimensionality, small sample sizes, and highly significant variables (1,2). The “glmnet” package provided a series of models, with the value of the tuning parameter λ inversely related to the model’s complexity and deviation. The number of nonzero coefficients dropped as the value of the invisible λ rose from left to right. Ten-fold cross-validation was performed to discover the ideal λ values, and a value lambda = 0.0266 with log (λ) = -3.626844 was selected using minimal criterion. However, the value’s findings may vary significantly depending on the time of analysis. As a result, 10-fold cross-validation was performed up to 100 times, and cross-validated errors were averaged. Next, the cuprotosis-related prognostic genes derived from LASSO regression analysis were used in multivariate Cox regression analysis. The signatures were created by combining several cuprotosis-related prognostic genes, followed by the determination of the Akaike information criterion (AIC) value for each individual cuprotosis-related prognostic gene. Following that, the best prognostic signature was developed based on the lowest AIC value with the best goodness of fit.

A total of 452 patients in the training set were split into low-risk (CRG score< median value) and high-risk (CRG score > median value) groups and then conducted Kaplan- Meier survival analysis based on the median risk score. The “ggplot2” R software was then used to conduct principal component analysis (PCA). Furthermore, the testing and all sets were classified into low- and high-risk groups, with each subgroup conducted to Kaplan-Meier survival analysis and the creation of receiver operating characteristic (ROC) curves.

Chi-square tests were employed to investigate the correlations between the CRG score and clinical variables (age, gender, tumor site, and TNM stage). We ran univariate and multivariate analyses on the training and testing sets to see whether risk scores were independent of other known clinicopathological variables. Furthermore, we conducted a stratified analysis to see whether the CRG score preserved its predictive potential in various subgroups based on age, gender, T stage, N stage, M stage, tumor stage, and tumor site.

CIBERSORT was utilized to measure the quantity of 22 infiltrating immune cells in diverse samples from low- and high-risk cohorts to assess the ratios of TIICs in the TME. We investigated the relationships between the CRG score and the percentages of 22 invading immune cells. We also utilized boxplots to compare the amounts of immunological checkpoint expression between the low- and high-score cohorts. We also looked at the links between the two risk categories and MSI.

The mutation annotation format (MAF) from the TCGA database was developed using the “maftools” R package to differentiate the somatic mutations of SKCM patients into high- and low-risk categories. We also computed the tumor mutation burden (TMB) score for each SKCM subject in the two classes. We estimated the half maximal inhibitory concentration (IC50) values of chemotherapeutic medications routinely used to treat melanoma using the “pRRophetic” software to investigate variations in the therapeutic effects of chemotherapeutic agents in patients in the two classes.

Based on the results of the independent prognosis study, the clinical parameters and risk score were utilized to create a prediction nomogram using the “rms” package. Each variable in the nomogram scoring system was assigned a score, and the final score was calculated by summing the scores from all variables in each sample (11). The nomogram was evaluated using time-dependent ROC curves for 1-, 3-, and 5-year survivals. The nomogram calibration plots were utilized to show the predictive value between the anticipated 1-, 3-, and 5-year survival events and the practical actual results.

We investigated the relationship between tumor mutational burden (TMB) and microsatellite instability (MSI) in 33 cancer species and determined whether they were upregulated or downregulated. Furthermore, in pan-cancer, a correlation study between risk score and TME as well as stemness indices was done.

DMEM medium supplemented containing 10% fetal bovine serum (FBS) was used to maintain A375 and HaCaT cell lines. These cells were cultured in a humidified chamber at 37°C with 5% CO₂. The extraction of total RNA from cells was done by using a Trizol reagent (Sangon Biotech, China) according to the instructions supplied by the manufacturer. Ultraviolet absorption spectrometry was used to assess RNA quantity and quality. Evo M-MLV RT Mix Kit (Accurate Biotechnology, Hunan, China) from total RNA (1 μg) was used for cDNA synthesis, Evo M-MLV RT Mix Kit (Accurate Biotechnology, Hunan, China) from total RNA (1 μg) was used. Quantitative PCR (qPCR) was performed using the SYBR Green Premix Pro Taq HS qPCR Kit (AG11718, Accurate Biotechnology, Hunan, China). Primers shown in Table S3 were used for amplifying products. The ΔΔcycle threshold method was used for quantifying the relative gene expression. GADPH acted as an endogenous control to normalize the expression level of relative gene expression. All results were obtained from at least three independent experiments with a minimum of three replicates per condition.

The K-M technique was used. To evaluate statistical significance, the log-rank test is utilized. R Computing Environment v4.1.2 was used for further analytical experiments. To study the prognosis of the estimate models for one, three, and five years, the ROC curve and AUC were exhibited using the “SurvivalROC” algorithm.

The Figures 1 and S1 depict the methodological approach used in this research. This research includes a total of ten CRGs. The SKCM cohort had a rather high mutation frequency, according to a summary study of the prevalence of somatic mutations in these 10 CRGs (Figure 1A). Of the 85 SKCM samples, 70 (85.35 percent) contained mutations in the CRGs (Figure 2A). CDKN2A exhibited the greatest mutation frequency (55%), followed by MTFF1, DLAT, DLD, GLS, PDHA1, LIAS, and others.

Following that, we investigated somatic copy number alterations in these CRGs and discovered common copy number alterations in all ten CRGs. GLS, MTF1, LIAS, FDX1, DLAT, PDHA1, CDKN2A, and PDHB all had widespread copy number variation (CNV) rises, whereas DLD and LIPT1 had CNV diminished (Figure 2B).

The locations of the CNV alterations in the CRGs on their respective chromosomes are depicted in (Figure 2C). We then evaluated by comparing the mRNA expression levels of SCKM and normal tissues and discovered that most CRG expression levels were significantly associated with CNV alteration. CRGs with low and moderate CNV, such as PDHA1, DLAT, GLS, and LIAS, were overexpressed in CRC comparison to those in normal melanoma samples, whereas CRGs with high CNV, such as PDHB, were significantly increased in melanoma samples (Figure 2D), implying that CNV may regulate CRG mRNA expression. While CNV can explain many observed changes in CRG expression, it is not the only factor involved in expression level regulation (12). Other factors that can influence gene expression include DNA methylation and transcription factors (13, 14). Our findings revealed a significant difference in the genetic landscape and expression levels of CRGs between melanoma and control samples, revealing a latent role for CRGs in melanoma tumorigenesis.

Seven hundred and eighty-one patients from four eligible melanoma cohorts (TCGA-SKCM, GSE19234, GSE65904, GSE78220, and GSE133713) were included in our study for further analysis to fully understand the expression pattern of CRG involved in tumorigenesis. Table S1 contains detailed information on the 781 melanoma patients. The prognostic values of ten CRGs in patients with SKCM were revealed by univariate Cox regression and Kaplan-Meier analysis (Table S4), and p< 0.05 was chosen as the filtering threshold. Following that, we ran a multivariate Cox regression analysis on five prognostic CRGs, three of which (FDX1, LIAS, and MTF1) were identified as independent potential predictors (Table S5). The comprehensive landscape of CRG interactions, regulator connections, and their prognostic value in melanoma patients was shown in a Cuprotosis network colored red (Figure 3A and Table S6). To investigate the expression properties of CRGs in melanoma, we used a consensus clustering technique on patients with melanoma based on the expression levels of the 10 CRGs (Figure S2). Our findings indicated that k = 2 showed up to be an efficient selection for sorting the entire cohort into subgroups A (n = 400) and B (n = 381) (Figure 3B). PCA analysis revealed significant differences in the transcription profiles of Cuprotosis between the subgroups (Figure 3C). The Kaplan-Meier curves revealed that patients with subtype A had a longer OS than patients with subtype B (log-rank test, p< 0.05; Figure 3D). Moreover, evaluations of the clinical and pathological features of various melanoma subtypes showed significant differences in CRG expression and clinical pathological characteristics (Figure 3E). As can be seen in Figure 3E, cluster A was associated with lower tumor location (p< 0.05), lower Stage (p< 0.05), and lower OS (p< 0.05) compared to cluster B.

Subtype A was significantly enriched in immunological fully-activated pathways, including cell cycle, p53 signaling pathway, pyruvate metabolism, and valine leucine and isoleucine degradation (Figure 4A; Table S7), according to GSVA enrichment analysis.

To examine the influence of CRGs in the TME of melanoma, we used the CIBERSORT algorithm (Table S8) to evaluate the correlations between the two subtypes and 22 human immune cell subsets of each melanoma specimen. We found significant differences in the infiltration of most immune cells between the two subgroups (Figure 4B).

Plasma cells, B cells memory, T cells CD8, CD4 memory-activated T cells, activated NK cells, monocytes, M1 and M2 macrophages, resting dendritic cells, and Eosinophils were clearly higher in subtype A than in subtype B, whereas B cells naive and T cells CD4 naive, T cells follicular helper, NK cells activated, M0 macrophages, mast cells activated had significantly lower infiltration in subtype A compared to subtype B. Similar fashion, analysis of two important immune checkpoints revealed that CLTA-4 and PD-L1 were more expressed in subtype B (Figure 4C, D).

We also used the ESTIMATE package to calculate the TME score (stromal score, immune score, and estimate score) for the two subtypes. Higher stromal or immune scores in the TME represented higher relative contents of stromal cells or immunocytes in the TME, whereas estimate scores indicated stromal or immune score aggregation in the TME. The results showed that patients with subgroup A had higher TME scores (Figure 4E).

We distinguished 409 Cuprotosis subclass DEGs using the R package “EdgeR” and conducted functional enrichment analysis (Figures 5A, B; Table S9) to investigate the potential biological behavior of each Cuprotosis pattern. These Cuprotosis subcategory genes were significantly enriched in biological processes associated with immune function (Figure 5A). KEGG analysis revealed an enrichment of immune and cancer-related pathways (Figure 5B), implying that Cuprotosis is important in the immune regulation of the TME. We, therefore, used univariate Cox regression to determine the prognostic value of 350 subcategory genes and filtered out 227 genes associated with OS time (p< 0.05), which were used in the statistical evaluation (Table S10).

To validate this regulation mechanism further, a consensus-clustering algorithm was utilized to split patients into 3 genomic subclasses based on prognostic genes, namely, gene subclasses A-C (Figure S3). Kaplan-Meier curves revealed that patients with gene subtype B had the lowest OS, while patients with gene cluster C had a better OS (log-rank test, p< 0.001; Figure 5C). Furthermore, Cuprotosis gene subtype B patterns were linked to an advanced stage and a greater risk of OS (Figure 5D). The three Cuprotosis gene subgroups revealed substantial changes in CRG expression, which was consistent with the Cuprotosis patterns (Figure 5E).

The CRG score was calculated using the DEGs. Figure 6A depicts the subject distribution in the two Cuprotosis subgroups, three gene subgroups, and two CRG score classes. First, we utilized R “caret package” to randomly assign subjects to train (n = 452) and testing (n = 191) groups in a 7: 3 ratio.

To determine the best prognostic signature, LASSO and multivariate Cox analyses were conducted on 227 Cuprotosis subgroup-related prognostic DEGs. Following LASSO regression analysis, 15 OS-associated genes were retained based on the least partial likelihood of deviance (Figures S4A, B).

We afterward conducted multivariate Cox regression analysis on 15 OS-associated genes using the Akaike information criterion (AIC) value to obtain six (SUSD2, SCYL1, KLF9, GSPT1, KRT73, ZNF780A), which included four high-risk genes (SUSD2, KLF9, KRT73, and ZNF780A) and two low-risk genes (SCYL1, and GSPT1) (Figure S4C). The CRG score was calculated using the findings of the multivariate Cox regression analysis as follows:

Risk score= 0.0290 ∗ SUSD2 − 0.0276 ∗ SCY L1 + 0.0262 ∗ KLF 9 - 0.0263 ∗ GSPT 1 + 0.0266 ∗ KRT 73 + 0.0265 ∗ ZNF 780A

The CRG score differed significantly amongst Cuprotosis gene subtypes. Genecluster A had the lowest CRG score, whereas genecluster C had the highest score (Figure 6B). More crucially, subgroup A had a much higher CRG score than subgroup B. Figure 6C depicts the patterns of risk scores in the two subgroups. Subjects with a CRG score less than the median risk score were classified as low-risk (n = 225), while those with a CRG score larger than the median risk score were classified as high-risk (n = 225). The distribution plot of the risk of CRG score demonstrated that as CRG scores rose, survival times fell while OS rates increased (Figure 6D, E). PCA analysis revealed distinct dimensions between the low- and high- CRG score groups (Figure 6F). The Kaplan-Meier survival curves demonstrated that patients with low scores had considerably better overall survival than those with high scores (log-rank test, p< 0.001; Figure 6G). Furthermore, the 1-, 3-, and 5-year CRG score survival rates were indicated by AUC values of 0.799, 0.717, and 0.710, respectively (Figure 6H). We computed CRG scores across the internal (testing set) and one external validation group (GSE65904) (Figures S5-6) to verify the prognostic effect of the CRG score.

The subjects were also divided into low- and high-risk groups based on the training set formula. The CRG scores, patient survival status, and PCA demonstrating the variation tendencies of the low- and high-risk groups are shown in Figures S5A, B, and S6A, B, respectively. Survival analysis demonstrated that the low-risk group had a considerably better prognosis than the high-risk group (log-rank; p< 0.001; Figures S5-6C). The CRG score still had reasonably high AUC values (Figures S5-6D) when the 1-, 3-, and 5-year prognostic prediction classification efficiencies were examined, demonstrating that the CRG score had an outstanding capacity to predict the survival of melanoma patients.

We investigated the link between the CRG score and several clinical features to determine the influence of the CRG score on clinical characteristics (age, sex, tumor location, and Stage). We found that patients in the stage IV subtype had substantially higher CRG scores than those in the stage I subgroup (p< 0.05; Figure S7A). To see whether this prognostic CRG score might predict OS independently in melanoma patients, we coupled clinical characteristics with the prognostic CRG score and ran univariate and multivariate analyses. The Stage and CRG scores in the training set indicated substantial differences, as shown in Figures S7B, C, with similar findings found in the testing (Figures S7D, E), and GSE65904(Figures S7F, G) groups.

Furthermore, a stratified analysis to determine whether the CRG score maintained its predictive ability in different subgroups, including age (= 60 and > 60 years), sex (female and male), tumor location (in-transit, local, primary, and regional), stage (stage I-II and stage III-IV), revealed that patients with high-risk scores had significantly lower OS compared to those with low-risk scores age (p< 0.01), sex (p = 0.0064 in women and p< 0.0001 in men), tumor location (p< 0.0001 in-transit, local, primary, and regional), stage (p< 0.001)(Figure S8).

We used the CIBERSORT method to examine the relationship between CRG score and immune cell abundance. The CRG score was favorably associated with M1 macrophages, plasma cells, activated memory CD4 + T cells, and CD8 + T cells, and negatively associated with nave B cells, M0 macrophages, follicular helper T cells, and gamma delta T cells (Figure 7A).

A high CRG score was linked to a high stromal score, immunological score, and estimate score (Figure 7B). We also looked at the link between the six genes in the suggested model and the number of immune cells. Most immune cells were shown to be highly associated with the six genes (Figure 7C). In addition, we looked at the relationships between immunological checkpoints and our risk model. Figure 7D indicates that 32 immunological checkpoints, including PD-L1 and CTLA-4, were expressed differently in the two groups.

Correlation studies indicated that a low CRG score was strongly connected to MSI-L status, while a high CRG score was related to MSI-H status (Figures 8A, B). We conducted survival analysis in the MSI-L and MSI-H groups to assess the impact of MSI status on OS in patients with melanoma. Although not statistically significant, the MSI group had a propensity for longer survival (p = 0.33; Figure S9).

Moreover, we also computed the tumor mutation burden (TMB) score for each SKCM subject in the two risk groups. Our study of the mutation data from the TCGA-SKCM cohort revealed a high TMB in the high-score group than in the low-score group (Figure 8C. The BRAF, PCLO, RP1, ANK3, APOB, DNAH7, ADGRV1, LRP1B, PTPRT, and THSD7B were the top 10 mutant genes in the high- and low-risk groups (Figures 8D, E). Patients with a low CRG score exhibited significantly greater rates of BRAF, PCLO, and LRP1B mutations than those with a high CRG score. The mutation levels of APOB and DNAH7, on the other hand, were exactly the reverse (Figures 8D, E). We next chose chemotherapeutic medicines, BRAF/MEK inhibitors (Dabrafenib, Trametinib) and immunotherapy drugs (pembrolizumab, nivolumab, ipilimumab) that are already used to treat melanoma to assess the sensitivities of individuals in the low- and high-risk categories to these medications. Interestingly, we found that the patients in the high CRG score group had higher IC50 values for Dabrafenib, pembrolizumab, nivolumab, ipilimumab, Bleomycin (Figure 9). However, patients in the high CRG score group had lower IC50 values for Trametinib.

Results for the other chemotherapeutic medicines are given in Figure S10. Together, these results showed that CRGs were related to drug sensitivity.

In this study, a nomogram integrating the CRG score and clinicopathological characteristics was developed to predict 1-, 3-, and 5-year OS rates in patients with SKCM (Figure 10). CRG score, age, and patient stage were among the predictors. Our AUC trials on the nomogram model revealed greater accuracy for OS at 1, 3, and 5 years in the training, testing, and one external validation set (Figures 10B–D). Our comparison of the nomogram’s prediction performance with that of the age, gender, tumor location, and stage in the three sets (Figure S11) revealed that the nomogram had 1-, 3-, and 5-year AUC values in the training set of 0.799, 0.717, and 0.710, respectively, whilst those of the age, gender, tumor location, and stage were lower (Figures S11A-C). The 1-, 3-, and 5-year AUC values of the nomogram in the testing set were 0.786, 0.722, and 0.657, respectively, whereas those of age, gender, tumor site, and stage were reduced (Figures S11D-F). Interestingly, the AUC values of the nomogram in two external validation sets (GSE65904) were greater than those of the age, gender, tumor site, and stage (Figures S11 G-I). The calibration graphs that followed showed that the suggested nomogram performed similarly in both the training and testing sets compared to an ideal model (Figures 10E–G).

We used pan-cancer analysis to assess the similarities and distinctions of the risk score model across malignancies (Figures 11A–E). TMB and MSI were carefully examined across malignancies. TMB was connected with the risk score in COAD, GBM, HNSC, KICH, LUAD, LUSC, READ, SKCM, STAD, and UCS (P<0.05, Figures 11A). MSI has a presence in LAML, LUAD, MESO, PRAD, SKCM, and STAD (P<0.05, Figures 11B). In addition, we computed the relationship between the risk score and 22 indices of immune cell infiltration and stemness. The outcomes are shown in Figures 11C–E.

RT-qPCR was used to detect the RNA level in human melanoma cells A375 and human normal cells HaCaT for validation work of the expression of RNAs identified in the risk model. A375 is a highly differentiated human malignant melanoma cell line and HaCaT is human normal skin immortalized keratinocyte-forming cell line from non-tumor tissues. The results of RT-qPCR in Figure 12 were consistent with biological analysis. We found that the expression of FDX1, LIAS, LIPT1, DLD, DLAT, MTF1, and CDKN2A were significantly lower in melanoma than in normal tissues. However, the expression of PDHA1, PDHB, and GLS was significantly higher in melanoma than in normal tissues.