Establishment and definition of the IRGP signature
A flowchart of the establishment and validation of 23 IRGPs is presented in Figure S1. The TCGA dataset was divided into a training dataset (n = 236) and a testing dataset (n = 235). 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 of > 0.5 were common among all datasets. After removing IRGPs with a score of 0 or 1 in less than 20% or more than 80% of the samples in the TCGA training and testing datasets, a total of 74,750 IRGPs had remained. Of these, 6,800 prognostic IRGPs were significantly associated with OS (FDR < 0.001), as determined with the log-rank test. Finally, 23 IRGPs were selected for the LASSO penalized Cox regression model (Figure 1a), including 46 unique IRGs, of which most coded for molecules related to antimicrobials and cytokines (Table 1). What’s more, we analyzed the Principal Component Analysis (PCA) of 46 unique IRGs by GEPIA web (gepia.cancer-pku.cn/), and found these genes expressed in TCGA tumor samples were independent parts, compared with normal skin that had sun explored or not (Figure S2).
The 23 IRGPs were used to calculate the risk score to predict the 5-year OS rate of each patient in the training cohort. 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 into the high or low immune risk group was -0.674 (Figure 2a). The data suggested that the high risk group had a higher risk index than the low risk group. A higher risk score indicates a higher number of deaths (Figure 1b, c), indicating that OS was significantly poorer for the high risk group than the low risk group (p < 0.001) (Figure 1e and Table S1). As shown in Figure 1d, the AUC values[20] of the 1-, 3-, and 5-year OS rates of the training cohort were 0.909, 0.901, and 0.912, respectively, and the C-index was 0.775 (95% confidence interval [CI] = 0.748–0.802) (Figure 4e). A nomogram of OS was created by combining all of the clinicopathologicalfactors such as age, sex, tumor stage and IRGPs risk score to predict the prognosis of CM (Figure 2c). As we can see, IRGPs risk score made the great contribution to the nomogram, and the 1-, 3-, and 5-year calibration curves (Figure 2b) demonstrated promising predictive ability of the nomogram.
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 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 (Figure 3a, b, 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), OS was expected to be better for younger patients (p = 0.029, Figure 3c). 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 3d). Finally, all patients were stratified by sex into a male dataset or a female dataset. As shown in Figure 3e, there was little difference in the OS rate between males and females (p = 0.543), possibly because of the small number of samples.
Collectively, the results showed that the predictive ability of the 23-IRGP signature was independent from other clinical parameters and was predictive of OS of CM patients.
Verification of the 23-IRGP signature to predict OS
To determine whether the 23-IRGP signature had consistent prognostic value in different risk groups, the TCGA test datasets GES65904 (n = 214), GSE59455 (n = 141), and GSE54467 (n = 79) were applied for external validation. The risk score of each patient was calculated with the same 23-IRGP prognostic signature. Then, based on the median risk score, the patients were assigned to the low or high immune risk group. Interestingly, OS was poorer for the high immune risk group (Figure 4a–d). The results of multivariate Cox regression analysis (Table 2) showed that the risk score of the 23-IRGP signature was an independent prognostic factor after adjustment for age, sex, and tumor stage in the testing dataset (hazard ratio [HR] = 2.08, 95% CI = 1.41–3.05, p = 0.0002), as well as the datasets GSE65904 (HR = 1.98, 95% CI = 1.49–2.62, p = 1.87E-06), GSE59455 (HR = 1.64, 95% CI = 0.94–2.82, p = 0.042), and GSE22153 (HR = 2.07, 95% CI = 1.16–3.69, p = 0.014). Finally, the C-index values for the testing, GSE65904, GSE59455, and GSE22153 datasets were 0.636 (95% CI = 0.585–0.687), 0.650 (95% CI = 0.609–0.691), 0.557 (95% CI = 0.508–0.606), and 0.610 (95% CI = 0.537–0.683), respectively (Figure 4e).
Comparison of the IRG prognostic signatures of CM
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 [21] and TCGA CM dataset were merged. Then, significantly different expressed immunity-related genes were selected for the next procedure. 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 5a). 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. The 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 (Figure 5b, c). Moreover, the high-risk group had a significantly poorer OS than the low risk group (p < 0.001) (Figure 5e). As shown in Figure 5d, 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 5E). Collectively, these results demonstrated that the prognostic signatures of the 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.
Immune cell infiltration,tumor micro-environment and 6 key genes analysis
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 within the different risk groups. A comparative summary of the CIBERSORT output resulting from these two risk groups is shown in Figure 6a. Immune cells, such as M0, M1, and M2 macrophages, plasma cells, activated CD4+ memory T cells, monocytes, CD8+ T cells, follicular helper T cells, and gamma delta T cells were enriched in different risk groups. As shown in Figure 6b and D, M0 macrophages (p = 0.004) and M2 macrophages (p = 0.003) were significantly highly expressed in the high risk group, while the abundances of M1 macrophages (p = 0.001), activated CD4+ memory T cells (p = 0.005), monocytes (p = 0.047), plasma cells (p = 0.011), CD8+ T cells (p = 0.028), follicular helper T cells (p = 0.017), and gamma delta T cells (p = 0.014) were significantly enriched in the low-risk group (Figure 6c, e–j). Then, we estimated tumor micro-environment (TME) in two different groups and found the high risk group had higher tumor purity with lower immune cells and stromal cells infiltration (Figure 7a). What’s more, we were seeing that 3 key targets of immunotherapy: PD1, PDL-1, CTLA-4 were highly expressed in low immune risk group (Figure 7b,c,d), which means the low-immube risk group may have better effect with immunotherapy. Kalaora S et al reported that that the PSMB8 and PSMB9 overexpression was predictive for better survival and improved response to immune-checkpoint inhibitors of melanoma patients[22], which were expressed highly in our study’s low immune risk group (Figure 7e,f). Interestingly, PRAME which as an independent biomarker of metastasis Uveal Melanoma[23] were also significant expressed in high risk group (Figure 7g). Previous results were the same in the Immunohistochemical results downloaded from The Human Protein Atlas dataset (THPA)(Figure 7n-w), as THPA has no result about CTLA4, we only present the other 5 genes expression in melanoma tissue. What’s more, in GEPIA, the patients with high PD1, PDL-1, CTLA4, PSMB8 and PSMB9 expression had better OS (Figure 7h-l). cBioPortal (https://www.cbioportal.org/) was used to investigate the mutation, according to the results, PD1, PDL-1, CTLA4, PSMB8 and PSMB9 had 5%, 1.9%, 1.6%, 5% and 4% probability to mutate(Figure S3a), respectively. Poorer OS occured if these genes with mutation (Figure S3b). in addition, in GEPIA analysis, PRAME had no significant effect in OS, which need further study to verify (Figure 7m).
Biological function analysis in high risk group divided by 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 with the use of genes that were more highly expressed in the high risk group than the low risk group. According to the results of GO analysis, these genes were positively correlated with skin-related biological functions, such as 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 of < 0.05 is presented in Figure 8. Information of every GO term is shown in Table S1. As shown in Figure S4, multiple melanoma progression-related pathways, including oxidative phosphorylation[24, 25], retinol metabolism[26-28], and ribosome[29, 30], were significantly upregulated in the high risk group (padj < 0.05). Collectively, these results provided evidence of molecular mechanisms affected by the IRGP signature and, thus, predictive of the prognosis of CM patients.