Figure 1 shows the whole analytical process of the study.

The GEO (https://www.ncbi.nlm.nih.gov/gds/) database provided the GSE22153, GSE22154, GSE46517, GSE54467, and GSE65904 datasets. GSE54467 was downloaded with log2-transformed and quantile-normalized matrices, while GSE22153, GSE22154, GSE46517, and GSE65904 datasets had to normalize, subsequently. In principle, genes with multiple probes were expressed on the basis of the median. The mRNA expression profiles were presented in the fragments per kilobase million (FPKM) format. The information used in this research satisfied the following criteria: (1) mRNAs with nonzero expression degrees occupied 75% of all samples; and (2) the patients had accurate follow‐up times. TCGA (https://portal.gdc.cancer.gov/) provided TCGA-SKCM dataset that was used as a training cohort and included mRNA expression profiles of 454 CM specimens and the related clinical follow‐up data.

GSE65904 datasets were used as a test cohort as these datasets contained DSS (210 patients) and DMFS (150 patients) information. Moreover, the GSE22153, GSE22154, GSE46517, and GSE54467 datasets were integrated as a verification cohort as these datasets contained single OS (236 patients) information. While integrating the GSE datasets, the batch effect was removed by the combat role of the R software package “sva”. Supplementary Table 1 shows the basic data of the datasets included in this study.

A total of 49 immune checkpoints were involved in this study, including the B7-CD28 family [CD274 (PD-L1), CTLA4, B7-H3, ICOSLG, TMIGD2, PD-L2, ICOS, PD-1, and HHLA2] (18, 19), the TNF superfamily (TNFRSF14, TNFRSF18, TNFRSF25, TNFRSF4, TNFRSF8, TNFRSF9, TNFSF14, TNFSF15, TNFSF18, TNFSF4, TNFSF9, CD40, BTLA, CD27, CD40LG, and CD70) (20), and some other immune checkpoint members (BTNL2, ADORA2A, CD160, CD200R1, CD200, CD244, CD44, CD28, CD80, CD48, CD86, FGL1, ENTPD1, HAVCR2, IDO2, KIR3DL1, IDO1, LAG3, LAIR1, NRP1, LGALS9, NCR3, and TIGIT) (21–23).

Applying the R package “psych”, “reshape2”, and “factoextra”, we performed PCA on four subsets of the validation dataset based on the expression of 34 immune checkpoint molecules. The results were visualized using R package “ggplot2”.

Based on the expression of immune-related genes (IRGs), patients were classified by unwatched cluster exploration. The consensus clustering algorithm (24) determined the number and stability of the clusters. These steps were carried out on the basis of the ConsensuClusterPlus package (25) and repeated 1,000 times to ensure the stability of the classification (maxK = 10, reps = 1000, pItem = 0.8, pFeature = 1, and clusterAlg = “pam”). The code for unsupervised cluster analysis and the input file are available as supplementary source files.

We analyzed the possible effect of immune checkpoint inhibition in patients with melanoma by Tumor Immune Dysfunction and Exclusion (26) (TIDE, http://tide.dfci.harvard.edu/). In the analysis, patients with TIDE scores greater than 0 were officially defined as nonresponders to immune checkpoint inhibitor therapy, while those with a TIDE score less than 0 were defined as respondents (26).

The difference of genes in cluster 1 and cluster 2 was analyzed by the edgeR package in R software (27). Modified p-values < 0.05 and |log2FC (fold variation)| > 1 (|log2FC > 1| and the modified FDR < 0.05) were set to identify significantly expressed RNAs. The ggplot2 package was used to plot a volcano map.

We carried out KEGG enrichment analysis and immune pathway enrichment analysis by using the ClueGO (28) plug-in in Cytoscape (version 3.8.2) and visualized it. In addition, significant GO or KEGG terms or genes were defined as those with the adjusted-p < 0.05 and at least three related mRNAs.

The immune score of every sample with R software were determined by the ESTIMATE algorithm, and the Wilcoxon test (29) was used to further compare the variations in the extent of immune cell infiltration among various cluster groups. The ratios of 22 and 10 immune cell subtypes were assessed by the CIBERSORT and MCPcounter package on the basis of the expression file, respectively (30). Further analyses adopted the samples with p < 0.05 in the CIBERSORT and MCPcounter exploration outcomes. The comparison of variations in immune cell subtypes among different cluster groups was made by the Mann–Whitney U test.

Univariate Cox analysis was employed to first calculate the prognostic value of every ICP by the R/survival package, and IRGs with p < 0.05 were chosen as the seed ICP for Cox LASSO regression. Then, prognostic signatures with the R packages “survminer”, “glment”, and “survival” were identified by multivariate Cox regression. The risk marks for every patient in the training group were calculated according to the following formula: risk mark = expGene1 × βGene1 + expGene2 × βGene2 + expGenen × βGenen (exp, prognostic gene expression level; β, multivariate Cox regression model regression coefficients) (31). All the samples were divided into high- and low-risk mark sets randomly with the mean risk mark as the cutoff value. The survival of every group was analyzed with the Kaplan–Meier method and log-rank test. ROC curves and the corresponding areas under the curve (AUCs) were produced with the R package “survivalROC”.

The ICPBS was validated in various clinical subgroups and histopathological subtypes. Whether the ICPBS is a separate risk factor for CM was evaluated by analyzing the univariate and multivariate Cox regression of the ICPBS and other clinicopathological elements. The indexes of the multivariate Cox regression model (through the R package “rms”, “Hmisc”, “lattice”, “Formula”, and “foreign”) were employed to build the nomogram model. Then, the total risk score was calculated on the basis of every predictor in the nomogram, and the patients were divided into two groups, with the middle risk mark as the cutoff point (32, 33). The C-index, calibration curve, and ROC analysis were adopted to determine the predictive precision and discriminative capacity of the nomogram. The test and validation cohorts adopted a similar analysis process.

The clinical effectiveness of the nomogram was determined with decision curve exploration through the quantification of the net benefits at different threshold probabilities in the cohort (34). The net benefit was calculated by subtracting the ratio of all patients with false-positive results from the ratio of patients with true-positive results and by weighing the relative harm of forgoing interventions by comparing with the negative outcomes of an unnecessary intervention (R package “rms” and “rmda”) (35).

The software GSEA v4.0.3 (www.broadinstitute.org/gsea) was used to perform gene set enrichment analysis (GSEA). All of the samples were divided into high- and low-risk score groups on the basis of the middle cutoff. We input the expression profiles of the mRNAs, groups of samples, and background files.

We used the website GEPIA2 (36) to analyze the variation in the expression of ICPBS molecules between CM tissues in TCGA and normal skin tissues in GTEx. We also downloaded two datasets of GSE15605 and GSE46517, which contain the expression matrix of CM tissue and normal skin tissue; then, we calculated the expression difference between the two groups by the Kruskal method and corrected the p-value by the Benjamini–Hochberg method (37).

All the analyses in this study were carried out through R software (version 3.6.3). The immune scores and the abundance of immune cell infiltration of the two groups were compared by the Wilcoxon test. Comparisons among three or more groups were made by one-way ANOVA and Kruskal–Wallis tests. A p‐value < 0.05 was regarded statistically significant. All the p-values in this study were corrected by the Benjamini–Hochberg method.

First, ICPs coincident in the TCGA, GSE22153, GSE22154, GSE46517, GSE54467, and GSE65904 databases were screened ( Figure 2A ). We merged the four datasets (GSE22153, GSE22154, GSE46517, and GSE54467) and removed the batch effect ( Figure 2B ). The PCA before and after the merging suggested that we fully removed the batch effect ( Figures 2C, D ). Based on the ICPs, TCGA was classified into two groups by unsupervised clustering ( Figures 3A, B ). According to the prognostic exploration of the two groups, cluster 1 displayed great benefits in survival and no recurrence ( Figures 3C, D ). A multivariate heat map containing 34 ICP expression levels and clinicopathological data is shown in Figure 3E . We found that 34 ICPs were significantly differentially expressed between the two clusters ( Figure 3F ).

Based on the results of TIDE analysis, there were significantly more responses to immune checkpoint inhibition therapy in patients in cluster 1 than in those in cluster 2 ( Figures 4A, B ). Besides, the KM analysis of patients who had received immunotherapy and chemotherapy showed that the response of patients in cluster 1 to immunotherapy and chemotherapy was significantly better than that of patients in cluster 2 ( Figures 4C, D ). Therefore, we named cluster 1 as immunotherapy-sensitive group (ISG) and cluster 2 as immunotherapy-tolerant group (ITG). We further explored the difference in the tumor immune microenvironment between the two groups. In general, the tumor microenvironment includes various cell types (immune, interstitial, and endothelial cells), inflammatory mediators, and extracellular matrix (ECM) molecules (29). For the hidden mechanism between the therapeutic response and tumor immune microenvironment, a set of analytic approaches associated with the immune profile were adopted. Among TCGA cohorts, the stromal, immune, and ESTIMATE scores of the ISG were significantly higher than those of the ITG ( Figure 4E ). According to the Kaplan–Meier exploration of the CM data in TCGA cohorts, OS was notably shorter in the low stromal, immune, and ESTIMATE mark groups ( Figures 4F–H ).

Since therapeutic response was closely associated with immune infiltration score, immune cell infiltration variations between ISG and ITG were analyzed. The proportions of 22 and 10 immune cell subtypes were assessed by the CIBERSORT and MCPcounter package on the basis of the expression file, respectively. For further analyses, we used the samples with p < 0.05 in the CIBERSORT and MCPcounter exploration outcomes. Most of the differential immune cells between ISG and ITG were found in the CIBERSORT analysis, which were predominantly M0 macrophages, M1 macrophages, follicular helper CD4 T cells, resting NK cells, resting dendritic cells, and CD8 T cells ( Figures 5A–C ). In the MCPcounter analysis, the proportions of CD8 T cells, B lineage, cytotoxic lymphocytes, monocytic lineage, and myeloid dendritic cells in the ISG were higher than those in the ITG ( Figures S1A–C ).

We analyzed the differences in 58,385 genes between ISG and ITG; we found that 1513 genes were upregulated in ISG, while 1,307 genes were upregulated in ITG ( Figures S2A, B ). Functional enrichment analysis was also carried out for the exploration of the GO database terms and KEGG pathways related to the 2,820 significantly differentially expressed genes. According to the outcomes, enrichment in immune response-related pathways included regulation of B-cell differentiation, T-cell receptor signaling pathway, antigen processing, and presentation of endogenous peptide antigen, neutrophil activation involved in immune response, and response to interferon-gamma ( Figure 6B ). Moreover, according to KEGG pathway functional enrichment, PD-L1 expression and PD-1 checkpoint pathway in cancer, T-cell receptor signaling pathway, antigen processing and presentation, Th1 and Th2 cell differentiation, and Th17 cell differentiation were the main pathways related to the genes ( Figure 6A ).

Univariate Cox survival analysis evaluated 34 ICPs in total, and 21 ICPs with p < 0.05 were screened and included in the follow-up analyses ( Figure S3A ). According to Figure S3B , seven ICPs (lambda value = 7) subject to multivariate Cox regression analysis ( Figure S3C ) were identified by LASSO regression analysis. Ultimately, we identified seven ICPs predicting CM patient prognosis, including CD40, CD80, CD86, ICOSLG, KIR3DL1, LAG3, and TNFSF4. The risk mark of the ICPBS was calculated using the following formula: risk mark = (−0.0017 × expression of CD40) − (0.7615 × expression of CD80) + (0.0092 × expression of CD86) − (1.0053 × expression of ICOSLG) − (0.9862 × expression of KIR3DL1) − (0.0081 × expression of LAG3) − (0.0450 × expression of TNFSF4). The patients in the training and the validation sets were divided into high-risk and low-risk groups, with the middle-risk mark as the cutoff point ( Figures S4A, B ). Moreover, in both the training set and the validation set, the expression of the seven ICPs in the low-risk group was significantly higher than that in the high-risk group ( Figures S4C, D ).

In the TCGA cohort, the patients in the high-risk group had poorer OS and RFS than those in the low-risk group (log‐rank test p‐value < 0.0001, Figures S5A, C ). Similar results were observed in the validation set ( Figure S5D ). Moreover, whether survival predictions made with the ICPBS were accurate in the TCGA and validation cohorts ( Figure S5B, E ) was determined by implementing ROC curve analysis to assess the AUC values for the 3‐year OS for TCGA cohort (AUC = 0.670) and the validation cohort (AUC = 0.784). In order to determine whether the ICPBS is equally effective for the prediction of DSS and DMFS in CM patients, a similar analysis was performed on the test cohort (n = 210). According to the Kaplan–Meier curves for DSS and DMFS, patients with high-risk marks had poorer DSS and DMFS than those with low-risk marks (log-rank test p-value < 0.001, Figures S5F, H ). Moreover, according to the outcomes, the AUC for 3-year DSS and DMFS was 0.674 and 0.722, respectively ( Figures S5G, I ).

As shown in Figure S6 , except for KIR3DL1, the other six ICPBS molecules were differentially expressed between CM and normal skin tissues. Specifically, LAG3, CD80, CD86, and TNFSF4 were upregulated ( Figures S6B, C, E, F ), while CD40 and ICOSLG were downregulated in CM tissues compared with normal skin tissues ( Figures S6A, D ). Surprisingly, similar results were observed in GSE15605 ( Figure S7A ) and GSE46517 ( Figure S7B ) datasets.

As the most important element, the pathological stage affects CM patient survival (38), and other factors, including sex and age, are also important (39); hence, CM patients in the training and the verification cohorts were stratified by three clinical features. In all subgroups, including stage I, stage II, stage III, and stage IV; men and women; and older (age ≥ 60 years) and younger (age < 60 years), the patients in the high-risk groups had shorter OS time than those in the low-risk groups ( Figure S8 , p < 0.05). Moreover, based on the results of TIDE analysis, we further showed that there were significantly more responses to immune checkpoint inhibitor therapy in patients with low risk than in those with high risk ( Figures S9A, B ). The KM analysis of patients who had received immunotherapy and chemotherapy in both groups showed that the response of the low-risk group to immunotherapy and chemotherapy was better than that of the high-risk group ( Figures S9C, D ).

The univariate and multivariate Cox regression models included the following clinical candidate predictors: age, sex, stage, and ICPBS ( Figures 7A, B ). ICPBS was an independent risk factor in patients with CM (HR = 3.517, p < 0.001). The indexes of the multivariate Cox regression model ( Figure 7C ) were adopted to build the nomogram model. Surprisingly, the C-index of the nomogram for predicting survival was 0.735 (95% CI, 0.733–0.737). According to the calibration curve, the predicted values conformed to the observed values of the probabilities of 3-year, 5-year, and 7-year OS ( Figure 7D ). In the end, the total risk marks were calculated on the basis of every predictor in the nomogram model. The AUC values for 3-, 5-, and 7-year survival using the predictive nomogram were 0.705, 0.711, and 0.702, respectively ( Figure 7G ). Besides, in the ROC analysis, the AUC value of ICPBS was higher than that of stage, and the AUC value of the model was higher than either the ICPBS and AJCC stage. The Kaplan–Meier analysis showed that patients with a high-risk mark had a poorer OS than those with a low-risk mark (log-rank test p-value < 0.0001, Figure 7E ). Finally, clinical decision analysis (DCA) showed that the clinical benefit rate of the model without stage alone was higher than that of the model without ICPBS alone. However, the clinical benefit rate of the complete model with these two factors was higher than that of the model without either of them ( Figure 7F ). These results suggest that ICPBS is superior to the traditional TNM staging system in predicting prognosis. The combined model, including both of these two factors, had the best prognostic ability and predictive stability.

To determine the robustness of this model, similar analyses were conducted in the validation sets. ICPBS was a risk factor in patients with CM (HR = 1.569, p = 0.007, Figures S10A, B ). Next, the indexes of the multivariate Cox regression model were adopted to build the nomogram model. The C-index of the nomogram for predicting survival was 0.757 (95% CI, 0.755–0.759). According to the calibration curve, the predicted values conformed to the observed values of the probabilities of 3‐year, 5‐year, and 7-year OS ( Figure S10C ). In addition, the AUC values shown by the total risk mark for 3-, 5‐, and 7-year survival were 0.810, 0.798, and 0.877, respectively ( Figure S10F ). Moreover, in the ROC analysis, the AUC value of ICPBS was also higher than that of stage, and the AUC value of the combined model was also higher than either of them. According to the Kaplan–Meier curves of OS, patients with high total marks had poorer OS than those with low total marks (log-rank test p-value < 0.0001, Figure S10D ). Ultimately, DCA also showed that the clinical benefit rate of the model without stage alone was higher than that of the model without ICPBS alone and that the clinical benefit rate of the combined model with these two factors was higher than that of the model without either of them ( Figure S10E ).

Twenty KEGG pathways related to the ICPBS, including the MAPK signaling pathway, natural killer cell-mediated cytotoxicity, the P53 signaling pathway, T-cell receptor signaling pathway, B-cell receptor signaling pathway, pathway in cancer, primary immunodeficiency, and intestinal immune network for IgA production ( Figure S11 ), were identified by GSEA.