The researchers were granted approval to conduct the research by Departmental Research Ethics Committee at the Xiangya Hospital, Changsha, China. The study protocol was approved by the institutional review board of Xiangya Hospital. All the procedures were performed in accordance with the Declaration of Helsinki and relevant policies in China.
Sources of melanoma patients
This was a retrospective study using public data for comprehensive analysis. In total, two independent data-sets from different high-throughput platforms were used in this study including TCGA (n=476), and Cabrita (GSE65904, n=214).  Gene expression data together with corresponding clinical information were accessed from Gene Expression Omnibus (GEO http://www.ncbi.nlm.nih.gov/geo/) and Firehose Broad GDAC portal (http://gdac.broadinstitute.org/, accessed on Mar 8, 2021). Totally, 684 samples were available for analysis.
Gene expression data processing
The original microarray datafile (GSE65904) was normalized using the robust multi-array average (RMA) method from the “affy” package. For each data set, the expression profile was converted from probe level to corresponding gene symbol based on each set of annotation files without further standardization; probe sets were selected on the basis of overall highest overall expression for each gene. For all data-sets, this study only selected patients with complete overall survival (OS) information.
Construction of prognostic immune related gene pairs (IRGPs) signature
The identification of prognostic IRGPs was performed as described previously.  The TCGA cohort were used for discovery and training the model. To construct immune related prognostic signature, we obtained 13,996 immune related genes (IRGs) from the ImmPort database (https://immport.niaid.nih.gov)  on 3/8/2021. A variety of immune-related genes were included including signal transduction, cell communication and immune response. IRGs was measured only on all platforms involved in this study and had relatively high variability (determined by median absolute deviation >0.5). Each IRGP was computed by pairwise comparison the gene expression level in a specific sample or profile. More specifically, in a pairwise comparison, the output is 1 if the first element is larger than the later one and 0 for the different order. After removing IRGPs with a small variation and imbalanced distribution (MAD =0), the remaining IRGPs were left and selected as initial candidate IRGPs for prognosis prediction. From the initial candidate IRGPs, an IRGP index (IRGPI) was constructed using Lasso Cox proportional hazards regression with 10-fold cross validation (glmnet package, version: 2.0–18) and 25 gene pairs were used to define the final model. To stratify patients into low or high-risk groups, the optimal cut off of the IRGPI was determined using time-dependent receiver operating characteristic (ROC) curve analysis (survivalROC, version 1.0.3) at 1 year in the training cohort for overall survival.
Validation of the IRGPs signature
The IRGPI model was further evaluated in the validation of melanoma patients by the log rank test. We then assessed the IRGPI based on other clinical factors in the univariate and multivariate cox proportional-hazards analysis.
Immune cells infiltration in bulk tumor GEPs
In order to study the enrichment of immune cells in different risk groups, we used CIBERSORT,  a machine learning method for predicting immune cell infiltration of bulk tumor GEPs. CIBERSORT used support vector regression to estimate the enrichment of various immune cell types in GEPs. For each sample, CIBERSORT quantified the relative abundance of 22 infiltrating immune cells, including T cells, B-cell macrophages and natural killer cells. Using Monte Carlo sampling, CIBERSORT calculates the deconvolution p value of each sample to provide the estimated confidence. The TCGA GEPs were uploaded to the online analytical platform CIBERSORT portal website (http://cibersort.stan ford.edu/), using the default signature matrix at 1,000 permutations.
Gene ontology (GO) and gene set enrichment analysis (GSEA)
GO analysis was performed for the prognostic immune signature using FunRich (http://www.funrich.org) . Gene set enrichment analysis  was conducted using the Bioconductor package “fgsea” with 10,000 permutations. The log2 fold change between the gene expression profiles of high-and low-risk groups was taken. Gene sets from high- vs low- risk groups were compared. Biological processes involved in this study were downloaded from the Molecular Signature Database (MSigDB C2 and C5 databases, version 7.2) (http://www.broadinstitute.org/gsea/msigdb/collections.jsp). FDR-adjusted P<0.05 or nominal (NOM) P<0.05 was used to select statistically significant gene sets.
All the statistical tests were performed using R (version 3.6.3, www.r-project.org). Student’s t-tests were used to compare the differences among groups. Cumulative survival time was calculated using the Kaplan-Meier method and the differences in survival curves were analyzed using the log-rank test from “survival” package (version: 2.42.3). Hazard ratios were calculated using the “survcomp” package  (version: 1.40.0). Univariate and multivariate analyses were conducted using the Cox proportional hazards regression model. For all tests, a pvalue of less than 0.05 indicated a statistically significant difference. Statistical significance is shown as *P<0.05, **P<0.01, ***P<0.001.