A Novel Signature of 23 Immunity-Related Gene Pairs Is Prognostic of Cutaneous Melanoma

In this study, we aimed to identify an immune-related signature for predicting prognosis in cutaneous melanoma (CM). Sample data from The Cancer Genome Atlas (TCGA; n = 460) were used to develop a prognostic signature with 23 immune-related gene pairs (23 IRGPs) for CM. Patients were divided into high- and low-risk groups using the TCGA and validation datasets GSE65904 (n = 214), GSE59455 (n = 141), and GSE22153 (n = 79). The ability of the 23-IRGP signature to predict CM was precise, with the stratified high-risk groups showing a poor prognosis, and it had a significant predictive power when used for immune microenvironment and biological analyses. We subsequently established a novel promising prognostic model in CM to determine the association between the immune microenvironment and CM patient results. This approach may be used to discover signatures in other diseases while avoiding the technical biases associated with other platforms.


INTRODUCTION
The incidence of cutaneous melanoma (CM) is increasing more rapidly than any other cancer, and accounts for about 132,000 new cases and 65,000 deaths worldwide annually (1). CM is the most lethal form of skin cancer and is a serious public health concern. The primary clinical feature of CM is early stage metastasis, which is one of the most significantly increasing cancers in the United States (2). Siegel et al. reported that in 2018, it had been 91,270 new cases and 9,320 deaths in the United States, owing to this disease (3).
As most diagnoses are made in the terminal stage, surgical results are often poor. At present, chemotherapy is the first line of treatment for CM (4), although many cases respond poorly to such regimens due to a high prevalence of adverse drug reactions and resistance to chemotherapeutic agents.
CM is associated with strong immunogenicity; thus, immunotherapy is a promising treatment alternative (5). Initial clinical trials using interferon-a (6) and high-dose interleukin-2 for advanced cases of CM (7) reported successful results. In addition, immune checkpoint inhibitors (ICIs), such monoclonal antibodies against cytotoxic T-lymphocyteassociated protein (CTLA-4) (8) and programmed cell death protein 1 (PD-1) (9), have provided meaningful results in the clinical outcomes against advanced melanoma, as demonstrated by improved survival and a greater curative effect for an increasing proportion of patients with CM.
However, despite broad efforts to identify novel prognostic biomarkers, the clinical behavior of CM remains unpredictable, rendering the prediction of survival time and tumor stage particularly difficult (10). Therefore, novel biomarkers and patient-tailored treatments are greatly needed, especially for patients at higher risk of recurrence. Although the immune system plays an essential role in the development and progression of CM (11), few immunity-related genes (IRGs) have been identified for use as tumor markers for prognosis. Nowadays, many recent CM investigations have several limitations such as studies only one or a few immunity-related biomarkers, small sample datasets, lack of follow-up validations or lack of detailed and comprehensive immunity-related studies. Moreover, many studies have reported that genetic mutations and the interactions between the tumor and its microenvironment can impact the biological features and malignant potential of CM. Considering that many immunityrelated treatments have shown significant therapeutic effects, identification of the interactions between cancer cells and the host immune response is especially important (12,13).
In this study, 23 IRGs associated with CM were identified from the whole transcriptome sequencing (RNA-seq) data retrieved from The Cancer Genome Atlas (TCGA) (14) and the ImmPort dataset (15). Then, three microarray datasets (GSE65904, GSE59455, and GSE22153) were selected from the Gene Expression Omnibus (GEO) database (16) to verify the usefulness of this 23 IRG pair (IRGP) signature for the prognosis of CM. Moreover, the possible relationships between clinical pathological factors and the prognostic signature were explored to validate the predictive efficacy and accuracy of the IRGP. Finally, analyses of immune cell infiltration, the tumor microenvironment, and biological functions of different risk groups based on the 23 IRGPs were performed.

CM Patient Data
In this retrospective study, four independent RNA-seq datasets and clinical data from diverse, high-throughput sequencing platforms were comprehensively analyzed. A CM dataset (n = 460) was downloaded from the TCGA data portal (https://portal. gdc.cancer.gov) and randomly assigned to a training dataset (n = 230) or a testing dataset (n = 230). In addition, the GSE65904 (n = 214), GSE59455 (n = 141), and GSE22153 (n = 79) datasets were downloaded from the GEO database (http://www.ncbi.nlm. nih.gov/geo/) for validation of the IRGP signature. In total, 905 samples were available for analysis. A file containing 1534 IRGs that was downloaded from the ImmPort database (https://immport.niaid.nih.gov) and the CIBERSORT analytical tool (https://cibersort.stanford.edu/) were used to determine an estimation of the abundances of 22 distinct cell types in a mixed cell population based on gene expression data. Immunohistochemical images of CM patients were downloaded from The Human Protein Atlas dataset (http:// www.proteinatlas.org/). All data was available for free from website.

Data Preprocessing
All data were preprocessed using R software (version 3.6.2). If more than one probe was matched to the same target gene, the average value of the probe was calculated to replace the expression level of the single gene. If there were multiple samples from the same patient, the average value of each gene was used as the expression level of that gene.

Establishment of Prognostic IRGPs Based on the TCGA Dataset
The TCGA CM dataset was randomly divided into a training group and a testing group by R package "caret" with the ratio of training group samples to test group samples is 0.5. The distribution of CM patients gender (p = 0.068), age (p = 0.047), clinical stage (p = 0.036), follow-up time (p < 0.0001), death rate (p < 0.0001) and the number of CM samples in the dichotomies was similar between the two groups (17,18). The GSE65904, GSE59455, and GSE2215 datasets were employed as the validation group. IRGs with relatively high variation (median absolute deviation >0.5) were extracted from the platforms, as described previously. For the pairwise comparison of a specific sample, two IRGs were paired off, and if the expression of the first IRG was higher than that of the second, the two IRGs were merged as an IRGP and assigned a score of 1; otherwise, a score of 0 was assigned. Then, IRGPs with score of 1 or 0 in over 80% of the specimens both in the training and testing groups were selected as potential prognostic IRGPs. Based on the results of a log-rank test, IRGPs with a false discovery rate (FDR) <0.001 (n = 23) were retained and entered into a least absolute shrinkage and selection operator (LASSO) penalized Cox regression model (iterations = 1000). The median value of the risk score was used as a cut-off to divide the patients into high-and low-risk groups. Next, a heat map, risk score map, state map of overall survival (OS), and 1-, 3-, and 5-year receiver operating characteristic (ROC) curves were created, and the concordance (C)-index was calculated. Then, the IRGPs were integrated with other clinical factors to construct a nomogram and a calibration curve. Finally, univariate and multivariate Cox regression analyses were performed to determine whether the 23 IRGPs were independent from other clinical parameters.

Verification and Assessment of the IRGP Signature for the Prediction of OS
The TCGA testing dataset and three microarray data files were selected to validate the signature comprised of 23 IRGPs. Every dataset was stratified into high-and low-risk groups based on the cut-off value of the prognostic signature. Next, the log-rank test and Cox analysis were performed, and a graph of OS was created to calculate the C-index between the high-and low-risk groups in each dataset. Finally, the signature of the 23 IRGPs was compared to the 1-, 3-, and 5-year area under the ROC curves (AUCs) and the C-indices.

Immune Cell Infiltration in CM
The CIBERSORT analytical tool (19) was used to explore the enrichment of immune cells in the two CM risk groups. CIBERSORT is a tool used for deconvolution of the expression matrix of immune cell subtypes based on the principle of linear support vector regression, which can estimate the enrichment of various immune cell types in CM. Based on the CIBERSORT results, a radar chart of significantly differentially expressed IRGs between the two risk groups was constructed. All procedures were performed using R software (version 3.6.2).

Statistical Analysis
All statistical analyses were performed using the software packages R (version 3.6.2; www.r-project.org) and Prism 8 (GraphPad Software Inc., San Diego, CA, USA). Training group and testing group were randomly divided by R package "caret". OS curves were plotted using the R packages "survival" and "survminer." A heat map of the IRGPs, risk score map, and OS status graph were created using the R package "pheatmap." A model of prognostic IRGPs was established using the R package "glmnet" (21). Univariate and multivariate Cox regression analyses were performed using the R package "survival." ROC curves were constructed using the R package "survivalROC." A nomogram and calibration curves were plotted using the R packages "rms," "foreign," and "survival." The C-index was computed using the R package "Hmsic." Immune cell infiltration was processed with the R packages "e701," "limma," and "fmsb" (22). The tumor environment plot, based on the R package "estimate" (23), and the expression levels of six single key genes were determine using the R package "ggpubr." GO and KEGG analyses were conducted using the R package "fgsea."

Establishment, Definition, and Genetic Variation of the IRGP Signature
A flowchart of the establishment and validation of the 23 IRGPs is presented in Figure 1. The TCGA dataset was divided into a training dataset (n = 230) and a testing dataset (n = 230) ( Table  S2). Filter analysis was applied to establish a prognostic model of 1,811 unique IRGs that were obtained from the ImmPort database (accessed on January 4, 2020). In total, 620 IRGs with a median absolute deviation >0.5 were common among the datasets. After removing any IRGPs with a score of 0 or 1 in <20% or >80% of the samples in the TCGA training and testing datasets, a total of 74,750 IRGPs remained. Of these, 6,800 prognostic IRGPs were significantly associated with OS (FDR < FIGURE 1 | A flowchart of the establishment and validation of 23 IRGPs. The TCGA data were divided into a training cohort (230), which was applied to identify potential IRGPs, and a testing cohort (230). A validate cohort of the datasets GSE65904, GSE59455, and GSE22153 was used to verify the 23 IRGPs, which were also compared to the IRGP-OBS and IRGs. 0.001), as determined by log-rank testing. Finally, 23 IRGPs consisted of 39 IRGs were selected for the LASSO penalized Cox regression model ( Figure 3A), including 39 unique IRGs, most of which encoded molecules related to antimicrobials and cytokines ( Table 1). Furthermore, we conducted Principal Component Analysis (PCA), bioligical function analysis and genetic and expression variation of the 39 unique IRGs using GEPIA web (gepia.cancer-pku.cn/), Metascape (https:// metascape.org/gp/index.html) and R package "RCircos" and "GenVisR". We found that the genes expressed in TCGA tumor samples were independent parts, compared to TCGA normal sample (only one sample), GTEx normal skin exposed samples or not to the sun (Figure 2A). What's more, Figure 2D showed 39 IRGs were positively correlated with cancer immunerelated biological functions, including cytokine-mediated signaling pathway, defense response to other organism, Influenza A, response to interferon-gamma, and type I interferon signaling pathway. We first summarized the incidence of copy number variations and somatic mutations of 39 IRGs in CM. The investigation of CNV alteration frequency showed that only 12 IRGs had alteration, and most were focused on the deletion in copy number, while CYBB, NGFR and BIRC5 had a widespread frequency of CNV amplification ( Figure 2B). The location of CNV alteration of 12 IRGs on chromosomes was shown in Figure 2C. Among 460 samples downloaded from TCGA-CM mutation dataset, 253 experienced mutations of 39 IRGs, with frequency 53.94%. We found that the LTBP2 exhibited the highest mutation frequency followed by PLXNB2, while almost antigen processing and presentation molecules (HSPA1A, ICAM1, and PSME1) as well as cytokines (CCL17 and LHB) did not show any mutations in CM samples ( Figure  2E). Further analyses revealed a significant mutation cooccurrence relationship between IKBKE and CYBB, IKBKE and CD8A, GPB2 and CD8A, GNLY and CD40, LYZ and TNFSF10, along with GPB2 and PSME1 ( Figure S1). To explore whether the above genetic variations influenced the expression of 39 IRGs in CM patients, we investigated the mRNA expression levels of genes between skin normal samples (GTEx skin normal sample and TCGA skin normal sample) and tumor samples, and found that the alterations of mutation could be the prominent impact factors resulting in perturbations on the 39 IRGs expression. Compared to normal skin tissues, 39 IRGs with high mutation obviously lower expression in CM samples (e.g., LTBP2, CD86, EDNRA, TRIM22, CYBB, STC1, GBP2, and RNASEL), and vice versa (e.g., LHB, PSME1, CCL17, ICAM1, ISG15, and BIRC5) ( Figure 2E and Figure S2). The above analyses presented the highly heterogeneity of genetic and expressional alteration landscape in 39 IRGs between normal and CM tissues, indicating that the expression imbalance of 39 IRGs played a important role in the CM occurrence and progression.
Twenty-three IRGPs were used to calculate a risk score to predict the 5-year OS rate of each patient in the training cohort. The analysis of the 5-year dependent ROC curve revealed that the best cut-off value of the 23 IRGPs to stratify patients in the training cohort and testing cohort into the high-or low-risk group was −0.674 ( Figure 4A). These data suggested that the high-risk group had a higher risk index than the low-risk group. A higher risk score means a higher number of deaths ( Figures 3B, C), indicating that OS was significantly poorer for the high-risk group than for the low-risk group (p < 0.001; Figure 3E). As shown in Figure 3D, the AUC values (24) for the 1-, 3-, and 5-year OS rates of the training cohort were 0.909, 0.901, and 0.912 (Table S3), respectively, and the C-index was 0.775 [95% confidence interval (CI) = 0.748-0.802] ( Figure 6E). Moreover, the AUC values for the 1-, 3-, 5-year OS rates of the testing cohort and TCGA dataset was also shown in Table S3. A nomogram of OS was created by combining all of the clinicopathological factors, including age, sex, tumor stage, and the IRGP risk score, to predict the prognosis of CM ( Figure 4C). The IRGP risk score made a major contribution to the nomogram, and the 1-, 3-, and 5-year calibration curves ( Figure 4B) demonstrated the promising predictive ability of the nomogram, moreover, the nomograms of TCGA-test dataset and TCGA dataset were shown in Figure S7. The 1-, 3-, and 5-year calibration curves of TCGA-test dataset and TCGA dataset were shown in Figures S9 and S10.
Univariate and multivariate Cox proportional hazards regression analyses of the TCGA dataset were performed to further assess the prognostic accuracy of the IRGPs for other clinical elements. The results of the univariate and multivariate Cox analyses showed that the signature of the 23 IRGPs and clinical factors, such as tumor stage, were indeed predictive of prognosis. However, although the IRGP signature was highly predictive of prognosis, the p-value was notably low (Figures 5A, B and Table 2).
Stratified analyses of patient age, tumor stage, and sex were also conducted. First, all patients in the TCGA training dataset were stratified by age into a young dataset (age ≤ 65 years) or an old dataset (age > 65 years), where OS was expected to be better for the younger patients (p = 0.029, Figure 5C). Then, all patients from the TCGA training dataset were further divided into an early onset dataset (stages I and II) or a later onset dataset (stages III and IV). Similar results were observed for the late dataset (p = 0.002, Figure 5D). Finally, all patients were stratified by sex into a male dataset or a female dataset. As shown in Figure 5E, there was little difference in the OS rate between males and females (p = 0.543), possibly due to the small number of samples.
Collectively, the results indicate that the predictive ability of the 23-IRGP signature was independent of other clinical parameters and was predictive of OS of CM patients.

Verification of Ability of the 23-IRGP Signature to Predict OS
To determine whether the 23-IRGP signature had consistent prognostic value in different risk groups, the validation datasets GES65904 (n = 214), GSE59455 (n = 141), and GSE54467 (n = 79) were applied for external validation. The risk score of each patient was calculated using the same 23-IRGP prognostic signature. Then, based on the median risk score, the patients were assigned to the low-or high-risk group. Interestingly, OS was poorer in the high-risk group ( Figures 6A-D). The results of multivariate Cox regression analysis ( Table 2) showed that, after    Figure 6E), the AUC values for the 1-, 3-, 5-year OS rates of the these datasets were also shown in Table S3. Moreover, the nomograms and the 1-, 3-, and 5-year calibration curves of three GSE validation datasets were shown in Figures S8, and S11-S13.
Immune Cell Infiltration, the Tumor Microenvironment, Potential of 23-IRGP as an Indicator of Response to Immunotherapy, and Analysis of Six Key Genes Reportedly, the infiltration of immune cells is associated with the prognosis of CM patients. The CIBERSORT analytical tool can be used to estimate the abundances of immune cell subsets and has been used in many previous studies of the cancer microenvironment. Therefore, the CIBERSORT analytical tool was applied to estimate the relative abundances of 22 different immune cells for each patient. A comparative summary of the CIBERSORT output resulting from the two risk groups is shown in Figure 7A.  Figures 7C, E-J). Then, we estimated the tumor microenvironment (TME) in the two groups and found that the high-risk group had a higher tumor purity with lower immune cells and stromal cell infiltration ( Figure 8A). Furthermore, as 23-IRGP had a potential of indicator of response to CM immunotherapy, the relationship between the 23-IRGP and ICIs, namely PD-1, PD-L1, and CTLA-4, were investigated. As shown in the Figure S14, the 23-IRGP was markedly negatively related with PD-1 and PD-L1 (rho = 0.321 and p < 0.001 for PD-1, and rho = 0.203 and p < 0.001 for PD-L1) ( Figures S14A, B), and positively correlated with CTLA-4 (rho = 0.085 and p = 0.145, without statistical significance). Moreover, three ICIs were found to be highly expressed in the low-risk group of 23-IRGP prognosis signature ( Figures 8B-D), and this result indicated that patients with low risk presented obviously higher expression levels of immune checkpoint genes (PD-1, PD-L1) compared with those in the high risk group (p < 0.001 for PD-1, and p < 0.001 for PD-L1) ( Figures 8B, C), which demonstrated that the PD-1 and PD-L1 are involved in better immunotherapy efficacy, and their high expression is associated with better prognosis. The effect of cross-talk between 23-IRGP and ICIs on CM patients' survival was shown in Figures S14D-F. According to the Meng Zhou et al. (25), we divided TCGA-CM patients into four clusters according to the connection of 23-IRGP and ICIs, and survival comparisons of three ICIs were presented among the four clusters. In PD-1 and PD-L1, Survival rate results showed that the 23-IRGP could significantly differentiate the result of patients with the same or similar levels of PD-1 and PD-L1 (P < 0.001, logrank test) (Figures S14D, E). Relative to other three clusters, patients who had low 23-IRGP value with high level of PD-1 or PD-L1 were likely to have best survival results. However, patients who had high 23-IRGP value with low level of PD-1 or PD-L1 expression tended to the poorest consequence compared with the other three clusters (Figures S14D, E). Meanwhile, patients who had low level PD-1 or PD-L1 with low 23-IRGP value had better survival outcomes than the patients that had low PD-1 or PD-L1 with high 23-IRGP value. Furthermore, no obvious statistical  significance was identified between the expression level of CTLA-4 and survival results for patients with 23-IRGP (P = 0.063, log-rank test) ( Figure S14F). Collectively, these investigations between the 23-IRGP and ICIs made us to speculate that the 23-IRGP may have a predictive ability of the response to CM immunotherapy. Kalaora et al. reported that that the among melanoma patients, overexpression of PSMB8 and PSMB9 was predictive for better survival and improved response to immune-checkpoint inhibitors (26), and these genes were highly expressed in the low-risk group ( Figures 8E, F). Interestingly, the PRAME gene, which is an independent biomarker of uveal melanoma metastasis (27), was also significant expressed in the high-risk group ( Figure 8G). These results correlated with the immunohistochemical results downloaded from The Human Protein Atlas dataset (THPA) ( Figures 8H-Q), which showed no results for CTLA-4, while the other five genes were expressed in melanoma tissue. Furthermore, in GEPIA, patients with high PD-1, PD-L1, CTLA-4, PSMB8, and PSMB9 expression showed better OS ( Figures S3A-E). cBioPortal (https://www.cbioportal.org/) was used to determine the mutation rates of the different genes. According to the results, the probability of mutation for PD-1, PD-L1, CTLA-4, PSMB8, and PSMB9 were 5%, 1.9%, 1.6%, 5%, and 4% ( Figure S4A), respectively. Poor OS was observed in cases where these genes were mutated ( Figure S4B). In addition, in the GEPIA analysis, PRAME had no significant effect on OS ( Figure  S3F). However, further study will be needed to verify this result.

Biological Function Analysis in the High-Risk Group Stratified by the 23-IRGP Signature
First, GSEA was used to investigate the results of the GO and KEGG pathway analyses between the high-and low-risk groups using genes that were more highly expressed in the high-risk group than the low-risk group. According to the GO analysis results, these genes were positively correlated with skin-related biological functions, including keratinization, epidermal cell differentiation, keratin filament, intermediate filament cytoskeleton, and skin development (padj < 0.05). A bubble graph of the 16 GO terms enriched in the high-risk group with a padj value < 0.05 is presented in Figure S5. Information on every GO term is provided in Table S1. As shown in Figure S6, several melanoma progression-related pathways, including oxidative phosphorylation (28,29), retinol metabolism (30)(31)(32), and ribosome (33,34), were significantly upregulated in the highrisk group (padj < 0.05). Collectively, the results obtained using the IRGP signature provide evidence of the molecular mechanisms affected by CM and, thus, the predictive power of this signature for the prognosis of CM patients.  could be defined as the time interval from TCGA sampling to patient death or last follow-up (35). As is shown in Figures 9A, B, IRGP-OBS prognostic signature divided the TCGA training dataset and testing dataset into a low-or high-risk OBS group, respectively, with cut off value −1.433 ( Figure 9D). Both in training and testing groups, high-risk group had a poor OBS than low-risk group. As shown in Figure 9C, the AUC values for the 1-, 3-, and 5-year OS rates of the training cohort were 0.946, 0.928, and 0.957, respectively, and the C-index was 0.751 (95% CI = 0.716-0.786) ( Figure 6E). It is worth noting that, group's immune cells infiltration of IRGP-OBS had a similar trend with 23-IRPG's ( Figure 9E). In this study, both IRGP prognostic and 23-IRGP prognostic models had predictive power of immune cells infiltration, with high AUC value and C-index value. The 23-IRGP prognostic signature was also compared with the prognostic signatures of individual IRGs. First, as the TCGA CM data had only one normal sample, the samples from the Genotype-Tissue Expression dataset (36) and TCGA CM dataset were merged. Then, significant differentially expressed IRGs were selected. Next, the LASSO penalized Cox regression model was applied to the TCGA clinical dataset, and 24 prognostic IRGs were selected for the final risk scoring model ( Figure 10A). Most of the 24 prognostic IRGs coded for molecules related to antimicrobials and cytokines. The IRGs significantly stratified the TCGA dataset patients into a low-or high-risk OS group. These data suggested that the high-risk group had a higher risk index than the low-risk group, as a higher risk score was associated with a higher risk of death (Figures 10B,  C). Moreover, the high-risk group had a significantly poorer OS than the low-risk group (p < 0.001) ( Figure 10E). As shown in Figure 10D, the AUC values of the 1-, 3-, and 5-year OS rates were 0.731, 0.760, and 0.749, respectively, and the C-index was 0.647 (95% CI = 0.612-0.682) ( Figure 6E). Collectively, these results demonstrate that the prognostic signatures of these IRGs had predictive ability but with a smaller AUC and lower C-index than the 23-IRGP signature, demonstrating that the 23-IRGP signature was the more precise predictive model in CM.

DISCUSSION
CM is a solid malignant tumor with strong immunogenicity with a rapidly increasing incidence rate worldwide. Since the approval of interferon-a for the treatment of CM in 1995, the potential of other immunotherapies have received much attention from researchers (37). Like with many other tumors, immune checkpoint (PD-1/PD-L1 and CTLA-4) blockade therapy has also become a target clinical. In 2011, ipilimumab, which targets CTLA-4, provided a major breakthrough in the clinical treatment of CM and was subsequently approved for marketing by the Food and Drug Administration (FDA). Ipilimumab is the first antibody drug to prolong the OS of patients with metastatic cancer (38). However, ipilimumab has been associated with some toxicity; thus, other immune surveillance sites have since been investigated, which has led to phase III clinical studies of the anti-PD-1 antibody drugs pembrolizumab and nivolumab. These drugs were approved for use by the FDA in September and December 2014, respectively. With lower anti-drug resistance and higher clinical safety, these anti-PD-1 antibody drugs offer hope to patients with advanced unresectable or metastatic CM (39,40). Given that the results of single antibody drugs are limited and the many links between the occurrence and development of CM, a multipleimmune therapy strategy may have more prospects (41). It is, therefore, necessary to develop a prognostic signature using IRGs.
Due to the technical bias of sequencing platforms, gene expression data must be preprocessed for standardization, which is particularly significant when establishing prognostic signatures. To achieve robust prognosis prediction without the technical bias associated with different platforms, prognostic signatures of IRGPs can be established by pairwise comparison, which does not require preprocessing for data standardization. IRGP scores are calculated based on the expression levels of IRGs in the same sample. As such, the | Tumor micro-environment (TME) and key genes expression in two different groups. In TME, "TumorPurity" is the percentage of tumor cells, "ImmuneScore" is the percentage of Immune cells, "StromalScore" is the percentage of stromal cells, "EstimateScore" is the percentage that merge the "ImmuneScore" and "StromalScore". as we can see, high risk group had higher tumor purity with lower immune cells and stromal cells infiltration (A). In low immune risk group, PD-1, PD-L1, CTLA-4, PSMB8, and PSMB9 were highly expressed (B-F). PRAME were significant expressed in high risk group (G). (H) IHC result of PD-1 protein in high risk group. Staining, not detected; intensity, negative; quantity, none; location, none. (I) IHC result of PD-1 protein in low risk group. Staining, low; intensity, weak; quantity, 75%-25%; location, cytoplasmic/membranous. (J) IHC result of PD-L1 protein in high risk group. Staining, not detected; intensity, negative; quantity, none; location, none. (K) IHC result of PD-L1 protein in low risk group. Staining, high; intensity, Strong; quantity, >75%; location, cytoplasmic/ membranous. (L) IHC result of PSMB8 protein in high risk group. Staining, low; intensity, moderate; quantity, <25%; location, cytoplasmic/membranous nuclear. (M) IHC result of PSMB8 protein in low risk group. Staining, high; intensity, Strong; quantity, >75%; location, cytoplasmic/membranous nuclear. (N) IHC result of PSMB9 protein in high risk group. Staining, not detected; intensity, negative; quantity, none; location, none. (O) IHC result of PSMB9 protein in low risk group. Staining, medium; intensity, moderate; quantity, >75%; location, cytoplasmic/membranous nuclear. (P) IHC result of PRAME protein in high risk group. Staining, medium; intensity, moderate; quantity, 75%-25%; location, cytoplasmic/membranous. (Q) IHC result of PRAME protein in low risk group. Staining, not detected; intensity, weak; quantity, <25%; location, cytoplasmic/membranous. ***p < 0.01 (t-test). prognostic signature is not only able to overcome the batch effects of different platforms, but also does not require the scaling and normalization of data. In CM prognostic model, based on the AUC and C-index values, the prediction capability of IRGPs is more promising when compared with prognostic checkpoint (AUC = 0.729) (42), prognostic DNA methylation (AUC = 0.822) (43), prognostic IRGs that require preprocessing for data standardization. Moreover, this approach has been reported to be robust in other cancer-related studies (44,45).
Given that the results of TCGA-CM dataset is made up of primary and metastatic samples, and the metastatic samples for sequencing were always from follow-up patients instead of initially diagnosed samples, IRGP-OBS prognostic model were established to make a predict comparison with 23-IRGP. As we expected, both IRGP-OBS model and 23-IRGP model had precious prediction capability with high AUC and C-index values. In addition, many researchers of CM also regarded OS as golden standard to evaluate the model predict power (43,46), and there is no significantly different trend in immune cell infiltration in our study. Collectively, IRGP model could be well applied to survival probability and immune cell infiltration of CM patients, providing reference for immunotherapy.
In the present study, an IRGP signature was established using a LASSO penalized Cox regression model to predict OS in CM patients. The prognostic signature of the 23 IRGPs consisted of 46 unique IRGs. Most of the genes in the immune signature encoded molecules related to antimicrobials and cytokines, which play important roles in the response to stimuli and the immune microenvironment. Many of these IRGs have be shown to be related to cancer development and prognosis, expression of serine/ threonine kinase 1 promotes melanoma metastasis (47), serum levels of C-C motif chemokine ligand 17 (CCL17) are an independent prognostic marker of distant metastasis of melanoma, and patients with 43% of patients with high CCL17 levels survived to 3 years (48). Singh et al. found that activation of intratumoral cluster of differentiation 40 induced T cell-mediated eradication of melanoma in the brain (49). Smith et al. discovered that endothelin 1 was enhanced in treated melanomas and conferred drug resistance via endothelin receptor type A (50). Ribonuclease L has been reported to interact with microRNA-146a as a sex-specific factor in melanoma (51), and semaphorin 7A has been found to reduce the pulmonary metastasis of melanoma (52). In this study, according to 23-IRGP signature, ICIs, including PD-1, PD-L1, and CTLA-4, were highly expressed in the low-risk group which had better survival. When exploring the relationship between PD-1/PD-L1 and the 23-IRGP value, the 23-IRGP signature presented significant relationship with PD-1/PD-L1 expression. Moreover, the interaction between ICIs and the 23-IRGP indicated a combined prognostic effect on patient survival. M2 macrophages have been shown to promote growth and are related to poorer OS in melanoma patients, while M1 macrophages support tumor destruction and antigen presentation (53). Yamaguchi et al. found that anti-PD-1 antibody (nivolumab) therapy increased the activated effector memory phenotypes of central memory T cells and subsets of CD4+ and CD8+ central memory T cells, as well as Th1 plus Thelper follicular 1 cells (54), indicating that these immune cells can prolong patient survival when they are activated. The results of the present study revealed a significant increase in the abundance of infiltrating M0 and M2 macrophages in the high-risk immune group, while the abundances of infiltrating M1 macrophages, activated memory CD4+ T cells, CD8+ T cells, follicular helper T cells, monocytes, plasma cells, and gamma delta T cells were greater in the low-risk immune group, which was found to have a better survival rate. In addition, both the mRNA analysis and immunohistochemistry results showed that the high-risk group had higher tumor purity and lower infiltration of immune cells and stromal cells. Meanwhile, low-risk group had higher immune According to the OS curve, OS was poorer for the high risk group as compared to the low risk group in the training cohort (p < 0.001). These results showed that the prognostic signature of the IRGs had good predictive power, but with a smaller AUC and lower C-index (0.610 in Figure 5E) than the 23 IRGPs signature. Therefore, the 23-IRGP was more precisely predictive for CM. checkpoint inhibitors indicating that patients in the low-risk group may have better outcomes with immunotherapy. Nowadays, ICIs provide a new way to cacer immunotherapy, when exploring the relationship 23-IRGP signature and ICIs, 23-IRGP signature showed closely connections with ICIs expression. Moreover, the relationship between 23-IRGP and ICIs indicated an integrated prognostic power on patient OS, which is associated with previous result that M1 macrophages, activated CD4+ memory T cells, monocytes, plasma cells, CD8+ T cells, follicular helper T cells, and gamma delta T cells infiltration and ICIs expression may make a difference on patients' OS and patients' immunotherapy effect. Collectively, it may suggested that the 23-IRGP may have a predictive ability of the response to CM immunotherapy (25).
In our previous study, overexpression of PSMB8 and PSMB9 was found to be predictive of better survival and improved response to immune-checkpoint inhibitors. This was reflected in the low-risk group in the present study, in which PSMB8 and PSMB9 were highly expressed, indicating that the low-risk group had better survival rate and immunotherapy effect. Moreover, this indicates that a mutation in these genes will result in poor OS. Unexpectedly, PRAME, which acts as an independent biomarker in uveal melanoma metastasis, was also significantly expressed in the high-risk group, indicating that PRAME may also be a biomarker in CM. However, further study will be needed to verify these results. Thus, our research outcomes were closely in line with those of previous studies, demonstrating the precise predictive power of our platform (55).
The 23-IRGP signature identified three pathways (i.e., oxidative phosphorylation, retinol metabolism, and ribosome) that were highly related to the invasiveness of melanoma, suggesting that a high-risk score was correlated with increased melanoma metastasis and poorer survival. These results indicated the capability of the IRGP signature for predicting tumor invasion in CM patients.
Nevertheless, there were some limitations to this study that should be addressed. First, the 23-IRGP prognostic signature was based on a retrospective study using the TCGA CM dataset and validated using three microarray datasets from the GEO dataset. Thus, these results should be validated against other datasets with different sample attributes in a prospective cohort. Second, as the 23 IRGPs were used to construct a prognostic signature model, different prognostic signature models are needed for comparison. Third, further validation of the 23 IRGPs by quantitative real-time polymerase chain reaction, western blotting, and immunohistochemical analyses will be needed before this approach can be applied clinically. Fourth, due to the fact that the relevant CM data (such as patients that received PD-1, PD-L1, and CTLA-4 treatment) cannot be obtained, the analysis of cross-talk between the signature and ICIs cannot be compensated systematically at present.

CONCLUSIONS
In conclusion, our findings indicate that our prognostic signature established using 23 IRGPs is a novel, promising model for predicting the prognosis of CM, indicating an association between the immune microenvironment and CM. This approach can be used to discover signatures in other diseases without the technical bias associated with different platforms.

AUTHOR CONTRIBUTIONS
Ya-NX conceptualized the project, all data analysis and wrote the first draft of the manuscript. Yi-NX, Z-CW, Y-ZM, and P-YW contributed to processing, analysis, and interpretation of the data. W-QT contributed to guide the data analysis, and manuscript writing. All authors contributed to the article and approved the submitted version.