3.1. Landscape of genetic variation of IMRGs in SKCM
In our study, we identified 85 IMRGs from the published literature and investigated their role in SKCM (Table S1) [34]. First of all, we determined the prevalence of somatic mutations in 85 IMRGs in SKCM. A total of 255 (54.6%) of the 467 samples experienced genetic changes in iron metabolism-related genes, including missense mutation, frame-shift deletion mutation, frame-in deletion mutation, splicing site mutation, and transcription start site mutation. The highest mutation frequency was found in HEPHL1, followed by HEPH and TF (Fig. 1A). Figure 1B shows the locations of CNV alterations of iron metabolism-related genes on chromosomes. Further analysis of 85 IMRGs showed that CNV mutations were common. HJV, BMP6, CUL1, HFE, NEO1, IREB2, FLVCR1, and EIF2AK1 showed widespread CNV amplification. In contrast, FTMT, NDFIP1, TF, LTF, ERFE, EPB42, B2M, TFR2, SLC39A8, BDH2, SRI and SLC25A37 have widespread CNV deletions (Fig. 1C). Moreover, we obtained RNA sequences of 812 normal skin tissue samples from UCSC Xena and combined them with 471 tumor samples obtained from the TCGA database. Next, we assessed possible changes in IMRG expression in normal skin tissue and melanoma samples. The results showed that some IMRGs were significantly reduced in normal samples, such as GLRX3, TFR2, STEAP2, TF, while others were significantly increased in tumor samples, such as NUBP1, SLC40A1, ATP6AP1, FLVCR1 and so on (***P < 0.001, **P < 0.01, *P < 0.05) (Fig. 1D). According to these analyses, the genomic and transcriptome landscape of regulatory factors of iron metabolism between normal and SKCM samples were significantly different and correlated. Hence, the expression alterations and genetic variation of IMRGs were crucial for regulating the occurrence and development of SKCM.
3.2 Identification of iron metabolism patterns mediated by 85 IMRGs
The GEO dataset (GSE65904) and TCGA-SKCM dataset with available survival data and clinical annotations were included in the cohort. We used the Kaplan-Meier survival method to understand the possible link between IMRGs with the OS of SKCM patients. It found that iron metabolism-related genes, such as ABCB6, ATP6AP1, ACO1, B2M, BMP6, BOLA2, CYBRD1, FLVCR1, HEPHL1, SLC22A17, TF, and so on, were related to the prognosis of SKCM (P < 0.05) (Figure S1, S2, S3). Next, we performed an unsupervised consensus analysis of all samples to stratify samples with different iron metabolism patterns based on the expression of iron metabolism-related genes. The result of k = 3 seems to be more accurate, which could divide all samples into three groups with less correlation between groups(Figs. 2A-B, Figure S4). Thus, we identified three distinct pattern clusters, including 187 cases in Cluster A, 315 cases in Cluster B, and 166 cases in Cluster C. Additionally, when we used PCA to identify changes in gene expression between the above three molecular subclusters, we found that the levels of iron metabolism-related genes could actually separate SKCM patients into three different molecular subclusters (Fig. 2C). The Kaplan-Meier analysis was used to compare the OS of the three clusters of patients to determine whether there was a correlation between the clustering and the clinical outcome. The results indicated that patients in Cluster A and Cluster C showed outstanding survival advantages, while Cluster B had the worst prognosis (P < 0.001) (Fig. 2D). The heatmap of the IMRG cluster evidenced that iron metabolism-related genes were differently expressed in three clusters (Fig. 2E). Besides, we observed significant differences in the expression of IMRGs among iron metabolism patterns.
Besides, GSVA was applied to identify the potential biological regulatory pathways among the different IMRG clusters in SKCM. According to the results of the two-by-two comparison between Cluster A and Cluster B, we can conclude that Cluster A is mainly concentrated in the toll-like receptor signaling pathway, natural killer cell media cytotoxicity, and chemokine signaling pathway (Fig. 2F). Based on the comparison between Cluster B and Cluster C, we identified that the JAK-STAT signaling pathway, primary immunodeficiency, and the intestinal immune network for IgA production were all active in Cluster C (Fig. 2G). In addition, through single-sample gene set enrichment analysis (ssGSEA) of SKCM, we could detect differences in immune cell levels across the three IMRG clusters discussed above in the immune microenvironment. We compared the levels of 23 immune cells in the three IMRG clusters to assess the association between iron metabolism-related genes and immune infiltration characteristics (Fig. 2H). The results showed that activated B cells, activated CD4 cells eosinophils, activated CD8 cells eosinophils, activated dendritic cells, eosinophilia, monocytes, neutrophils, and plasmacytoid dendritic cells were significantly different between different iron metabolism patterns. These results indicated that these iron metabolism modification patterns had different characteristics of immune infiltration, which may guide immunotherapy in patients with SKCM.
3.3 Iron metabolism phenotype related DEGs in SKCM
Based on the expression of IMRGs, a consensus clustering algorithm classified SKCM patients into three iron metabolism modification patterns, the potential genetic alterations and expression perturbations within these patterns were not well known. Therefore, we further examined the potential changes in iron metabolism-related transcriptional expression in three iron metabolic modification patterns in SKCM. We made the criterion adjusted P-value < 0.001 to identify overlapping differentially expressed genes (DEGs). A total of 703 DEGs that represented the critical distinguishing index of the three IMRG patterns were considered to be signatures of iron metabolism and illustrated by Venn diagrams (Fig. 3A, Table S3). Next, KEGG pathway enrichment and GO functional enrichment analysis were performed on the above DEGs. KEGG enrichment analysis showed that such DEGs mainly act in cytokine − cytokine receptor interaction, chemokine signaling pathways, NF − kappa B signaling pathway, and B cell receptor signaling pathway (Fig. 3B-C). GO enrichment analysis of these signature genes was mainly riched in T cell activation, leukocyte cell − cell adhesion, regulation of T cell activation, lymphocyte differentiation, and regulation of leukocyte cell − cell adhesion (Fig. 3D-E). These results provide further evidence that overlapping DEGs were characterized by iron metabolism modification and immune microenvironment regulation, which could be regarded as iron metabolism characteristics. In our unsupervised consistent cluster analysis, three transcriptome phenotypes were identified as gene clusters A, B, and C based on these 703 most representative iron metabolism patterns-related signature genes (Fig. 3F, Figure S5). Furthermore, survival analysis of SKCM samples revealed significant prognostic differences among three iron metabolism-related gene signatures. It has been found that gene cluster B is associated with better survival outcomes, whereas gene clusters A and C are associated with poorer outcomes (Fig. 3G). According to heatmaps based on different subclusters of DEGs and different clinical characteristics, the highest expression of DGEs was in cluster B, followed by cluster A, and then cluster C (Fig. 3H). Further analysis was conducted to visualize the expression of iron metabolism-related genes between different gene clusters using a boxplot, and we observed that iron metabolism-related genes were differentially expressed between the three gene clusters (***P < 0.001) (Fig. 3I).
3.4 Construction of the IMRG scores
While we identified the role of iron metabolism in prognosis, and the modulation of immune infiltration in SKCM patients, these analyses were based on only the patient populations and can’t accurately predict patterns of iron metabolism in individual tumors. Therefore we devised a scoring system called iron metabolism-related genes score (IMRG Score), which quantifies iron metabolism patterns in individual SKCM patients based on the IMRGs, to quantify the iron metabolism modification pattern of each SKCM patients. A Sankey diagram was used to illustrate the workflow of iron metabolism score construction because of the complexity of its quantification. Hence we further demonstrated the distribution of SKCM patients between IMRG subclusters, gene subclusters, IMRG score groups, and survival status. 47% SKCM patients were in the IMRG cluster B, 40% SKCM patients were in the gene cluster A and 54% SKCM patients were in the high IMRG score group (Fig. 4A). In addition, we attempted to determine the ability of IMRG scores to predict survival outcomes by classifying patients into low-score or high-score subgroups. As expected, patients with low IMRG scores were significantly associated with better outcomes in the cohort (P < 0.001, Fig. 4B). From this, it is clear that there is a close association between a low IMRG score and a good prognosis for SKCM patients. Besides, there was a significant difference in the IMRG score between IMRG clusters and gene clusters, with subcluster B of IMRG clusters and subcluster C of gene clusters possessing high IMRG scores (Fig. 4C and D).
Next, we investigated whether the IMRG scores for OS was an independent prognostic factor. Based on the results of multivariate Cox regression analysis, it was found that the IMRG score was an independent prognostic factor as well as age in SKCM patients (Fig. 4E). Based on all independent factors, a nomogram was developed to predict the 1-, 3-, and 5-year OS. Figure 4F illustrates that IMRG scores contribute the most to the total score. As the total score increased, the 1-, 3-, and 5-year OS rates decreased. As a result, these findings strongly suggest that the IMRG score can describe the pattern of iron metabolism modification and predict the future prognosis for patients with SKCM.
3.5 Analysis of correlations between IMRG scores and different clinical features
In this study, we extended our investigation into the association between IMRG scores and different patient characteristics such as age and gender. Based on different clinical characteristics, we divided SKCM patients into different subgroups, such as patients > 65 years old and > 65 years old, males and females. Furthermore, we compared the overall survival between patients with different IMRG scores and found that the OS was better for patients with a low IMRG score than for those with a high score in the age and gender subgroups (Fig. 5A-D). Besides, SKCM patients were divided into two groups based on the optimal cutoff for TMB (Tumor mutation burden) by using the R package 'survminer'. The result shows that the optimal cutoff for TMB is 3.11 (Figure S6). Furthermore, the overall survival was assessed in each group. It was found that patients with a high TMB had a longer survival time than those with a low TMB (Fig. 5E). Moreover, a combined analysis of the TMB and IMRG scores found that patients with a low IMRG score and a high TMB score had a better OS rate (p < 0.001) (Fig. 5F). It seems reasonable to conclude that iron metabolism-related genes and TMB are somehow related and together influence prognosis in SKCM patients.
3.6 The role of IMRG score in predicting immunotherapeutic benefits
In recent years, immunotherapy has been widely used in skin cutaneous melanoma cancer[35]. Among them, immune checkpoint inhibitors have made rapid progress in SKCM. Programmed cell death protein 1 (PD-1), programmed cell death ligand − 1 (PD-L1), and cytotoxic T-lymphocyte antigen-4 (CTLA-4) are the three main practice pathways. Anti-CTLA4 combined with anti-PD-1 checkpoint inhibitors has improved 5-year relative survival in patients with metastatic melanoma[36]. However, not all patients benefit from this therapy, and patients with certain tumors can cause serious immune adverse events. We first analyzed the differences between ICB genes (PDCD1, CD274, CTLA4) in the high and low IMRG score groups. We observed that the expression level of PD-L1, PD-1, and CTLA4 were significantly increased in the low IMRG score group (Fig. 6A-C). Next, we used the Cancer Immunome Atlas Database (TCIA), which provides results of comprehensive immunogenomic analyses and immunophenotype score (IPS) of next-generation sequencing data (NGS) data for SKCM patients, to evaluate the correlation between IMRG score and the immunotherapeutic response of SKCM. Consequently, we found that IMRG score was significantly associated with prognosis in melanoma patients treated with anti-PD1, anti-CTLA4, and anti-PD-1 in combination with anti-CTLA4 immunotherapy (Fig. 6D-F). At the same time, an SKCM patient with a low IMRG score had better efficacy than those with a high IMRG score, especially in the anti-PD-1 combined with the anti-CTLA4 group (p < 2.22e − 16).
Then, we analyzed the relationship between IMRG score and OS in the external validation cohorts, to further evaluate the predictive effect of IMRG score. We combined the two independent cohorts (GSE22153 and GSE19234). As we expected, the OS was better for patients with a low IMRG score (Log-rank test, P = 0.012). Besides, we analyzed the fraction of clinical outcome of SKCM patients in low and high IMRG score groups and found that 82% SKCM patients in the high IMRG score group had a poorer survival, indicating that the low IMRG score group had a better prognosis than the high IMRG score group (P = 0.0013).
In conclusion, these findings indirectly indicated that the representation of tumor iron metabolism plays a central role in mediating the immune response. Meanwhile, our results suggested that IMRG scores were associated with immunotherapy response in SKCM patients and may serve as a further predictor of patient outcomes.