The flow chart of the present study is shown in Scheme 1. The fragments per kilobase of per million (FPKM) of CM transcriptome, lncRNA counts data and corresponding clinical data of CM were downloaded from TCGA database (https://tcga-data.ncinih.gov/tcga/). The transcriptome samples included 471 tumor samples and 1 normal sample. Among them, 461 clinical samples were included for prognosis analysis because 10 patients had missing survival time. The clinical characteristics of the 461 patients are displayed in Supplementary Table S1. In addition, 233 normal skin samples with transcriptome data were downloaded from GTEx database (https://gtexportal.org/home/), then normalized and processed with TCGA biolinks package (22). The batch effects between TCGA and GTEx data were eliminated by R software package “limma”. The final datasets contained RNA expression profile of 471 CM samples and 234 normal samples. Sixty ferroptosis-related genes (Supplementary Table S2) that have been validated in previous studies (23–25) were obtained from the FerrDb database (26).

To identify the ferroptosis-related differentially expressed genes (FDEGs) between tumor and normal samples, we analyzed the datasets using the classical Bayesian algorithm in the “limma” package in R (27). The ggplot2 package was used to create heatmap plots (28). We performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses to evaluate the biological processes and pathways associated with these FDEGs. Analyses of Pearson correlation between lncRNAs and FDEGs were performed to identify ferroptosis-related lncRNAs. The association was considered significant if the square of the correlation coefficient |R²| > 0.4 and P < 0.001. Cytoscape software 3.8.0 was used to visualize the network of gene-lncRNA co-expression (29). All ferroptosis-related lncRNAs in the co-expression network were included in the differential expression analysis. Ferroptosis-related differentially expressed lncRNAs (FDELs) (false discovery rate [FDR] < 0.05 and log2FC ≥ 1) were identified using the “limma” package and illustrated in heatmap and volcano plots

We first performed univariate Cox regression to identify candidate FDELs that were significantly associated with overall survival (OS). We then performed least absolute shrinkage and selection operator (LASSO) regression to identify the key FDELs based on the optimal lambda value. We performed stepwise multivariate Cox regression including these key FDELs to construct a risk score. The risk score was calculated as follows: risk score =∑i=1nβi∗ (expression of lncRNA_i), where n is the number of key lncRNAs and β is the regression coefficient. Using the median risk score as a cutoff value, we classified all patients as high-risk (≥median risk score) or low-risk (<median risk score).

All 461 CM samples were randomly divided into a training cohort (n = 231) and a test cohort (n = 230) using the caret package in R. The detailed baseline characteristics of the training and testing cohort showed no statistical significance (Supplementary Table S3). The receiver operating characteristic (ROC) curves were applied and area under the curves (AUCs) were calculated to evaluate the accuracy of the prognostic lncRNA signature in the training, test, and entire cohorts. The risk level distribution and survival status distribution were illustrated. Kaplan-Meier analyses and log-rank tests were performed to compare survival differences between the high- and low-risk groups. In addition, principal component analysis (PCA) (30) and t-distributed stochastic neighbor embedding (t-SNE) (31) were used for dimensionality reduction and clustering visualization based on the risk score.

To determine whether the risk score was an independent predictor of survival, we performed univariate and multivariate Cox regression analyses, which included the risk score, patient age, patient sex, melanoma stage, and TNM stage (T: tumor size; N: lymph node involvement; M: metastasis), in the entire cohort. Hazard ratios (HRs) and 95% confidence intervals (CIs) for each predictor were calculated. AUCs for 3-year OS predicted with the independent risk factors as shown by multivariate Cox regression were calculated to evaluate the predictive value of the lncRNA signature. We also developed a nomogram that included the abovementioned risk factors for the practical prediction of 1-, 3-, and 5-year OS. The performance of the nomogram was evaluated using a calibration curve analysis and C-index.

We performed gene set enrichment analysis (GSEA) (32) to identify significantly differentially enriched pathways between the high- and low-risk groups. Gene expression was used as a phenotype label (33). The number of random sample permutations was set at 1,000. Statistical significance was set at a |normalized enrichment score| >1, nominal P-value <0.05, and FDR q-value <0.25 (34).

The tumor purity, ESTIMATE score, immune score, and stromal score for each melanoma sample were calculated using the ESTIMATE algorithm (17). In addition, the single-sample GSEA (ssGSEA) (35), TIMER, CIBERSORT (36), MCPcounter, QUANTISEQ, XCELL, and EPIC algorithms were used to compare cellular components and immune responses between high- and low-risk groups. The differences in immune cell infiltration were shown in heatmap, violin, and radar plots. Furthermore, the ssGSEA scores of the 16 most commonly involved immune cells and 13 most commonly involved immune functions in TME were compared between the 2 groups.

To predict the efficacy of ICIs, we systematically searched the ICI-related gene expression profiles and identified 44 ICI-related genes that might be correlated with ICI response. Among these genes, PD1 and CTLA4 were previously reported to be crucial targets of ICIs in CM (6, 37). The association between the expression of the PD1 and CTLA4 genes and the risk score was analyzed using data from The Cancer Immunome Atlas data (https://tcia.at/) to investigate the role of ferroptosis in the efficacy of ICIs.

The somatic mutational profile of CM was downloaded from the TCGA. The quantity and quality of the gene mutations were analyzed in the 2 groups with the “Maftools” package in R (38). In addition, we used the “pRRophetic” package in R to evaluate the half maximal inhibitory concentration (IC50) of common chemotherapy and molecular drugs.

All statistical analyses were performed using R version 3.6.1 (Institute for Statistics and Mathematics, Vienna, Austria; https://www.r-project.org) and its (packages: limma, impute, bioconductor, ggplot2, rms, glmnet, preprocess Core, forest plot, survminer, survival ROC, beeswarm, etc.) as described previously (39). Normally and non-normally distributed variables were analyzed using the unpaired Student' t-test and the Wilcoxon test, respectively. The Benjamin-i-Hochberg method was used to identify FDELs based on FDR. All statistical tests were 2-sided with P ≤ 0.05 being statistically significant.

We identified 53 FDEGs (25 upregulated and 28 downregulated) between the normal and tumor samples (Figure 1A). The results of GO (Figure 1B) and KEGG (Figure 1C) analyses indicated that FDEGs were mainly enriched in ferroptosis-related biological processes (e.g., response to oxidative stress) and signaling pathways (e.g., ferroptosis). Pearson correlation analysis identified 387 ferroptosis-related lncRNAs (326 upregulated and 61 downregulated) (Supplementary Table S4). The ferroptosis-related lncRNAs were visualized in a gene-lncRNA co-expression network (Figure 1D). As shown in the heatmap (Figure 1E) and volcano plot (Figure 1F), there were 161 FDELs (130 upregulated and 31 downregulated).

The univariate Cox regression analysis revealed that 50 of the 161 FDELs were significantly associated with OS (Figure 2A). LASSO regression revealed 8 of these 50 lncRNAs to be key to prognosis (Figures 2B,C). Stepwise multivariate Cox regression finally identified 5 lncRNAs (AATBC, AC145423.2, LINC01871, AC125807.2, and AC245041.1) for construction of the prognostic lncRNA signature (Figure 2D). Among these 5 lncRNAs, only LINC01871 was a protective factor, whereas the other 4 were risk factors (Figure 2E). The results of the univariate and multivariate Cox regression analyses of the 5 lncRNAs in the training cohort are shown in Supplementary Table S3. The risk score was calculated as follows: Risk score = 0.147 × expression of AATBC + 0.175 × expression of AC145423.2 - 0.091 × expression of LINC01871 + 0.165 × expression of AC125807.2 + 0.376 × expression of AC245041.1. All patients were divided into high- or low-risk group based on the median risk score of 0.995.

The AUC of the lncRNA signature for predicting 5-year OS in the TCGA training cohort was 0.743 (Figure 3A), indicating that the lncRNA signature offered accurate prediction. The mortality rate increased and survival time decreased as the risk score increased (Figures 3B,C). The OS rate of the high-risk group was significantly poorer than that of the low-risk group (Figure 3D). The ability of the 5-lncRNA signature to discriminate between high- and low-risk patients was validated in the testing cohort (Figures 3E–H) and entire TCGA cohort (Figures 3I–L). PCA and t-SNE analyses confirmed that the risk score had reliable clustering ability in the training cohort (Figures 3M,N), testing cohort (Figures 3O,P), and entire cohort (Figures 3Q,R).

The clinicopathological features of all 471 CM patients are shown in Supplementary Table S5. The HRs of the risk score in the univariate and multivariate Cox regression analyses were 1.499 (95% CI: 1.302–1.725, P < .001) and 1.338 (95% CI:1.154–1.552, P < .001), respectively (Figures 4A,B), indicating that the lncRNA signature was an independent predictor of OS.

We developed a nomogram that included the risk score and traditional clinicopathological factors for clinical use. The 1-, 3-, and 5-year survival probabilities for CM patients predicted with the nomogram were 92.7%, 57.2%, and 40.2%, respectively (Figure 4C). The AUC of the risk score for predicting 3-year OS (0.713) was greater than the AUCs of other clinicopathological features for predicting 3-year OS (684–0.693), indicating that the risk score had a better predictive ability (Figure 4D). In addition, the AUC of a prediction model that included the statistically significant variables in the multivariate analysis (i.e., risk score, age, T stage, and N stage) was 0.804, indicating that the novel risk score significantly increased the predictive accuracy of nomogram construction. The calibration curve analysis demonstrated excellent agreement between nomogram predictions and actual observations, indicating that the nomogram was stable and accurate in predicting the prognosis of CM (Figure 4E). Finally, the comparison of signature with other models showed that our AUC was superior to the existing signatures (Figure 4F). The calculation of C-index of Nomogram and signature were performed to reflect the validity of the nomogram and model, and two results are respectively 0.696 and 0.739 (Figure 4G).

The results of GSEA showed that ferroptosis, apoptosis, immune system process, immune response, antigen processing and presentation, T-cell receptor, B-cell receptor, and natural killer cell-mediated cytotoxicity pathways were inhibited in the high-risk group (Figures 5A–H). The lower level of ferroptosis and immune response in the high-risk group suggested that the association between ferroptosis and the immune TME play an important role in CM prognosis; therefore, we further compared immune cell infiltration and functions between the high- and low-risk groups.

Significant differences in immune cell infiltration between the high- and low-risk groups are illustrated in Figure 6A. Compared with the low-risk group, the high-risk group had significantly higher tumor purity (Figure 6B), but lower ESTIMATE scores (Figure 6C), stromal scores (Figure 6D), and immune score (Figure 6E), indicating that the high-risk group had a lower immune response. The high-risk group also had significantly fewer CD8+ T cells, CD4+ memory-activated T cells and M1 macrophages, but more M0 and M2 macrophages (Figure 6F). In addition, a comparison of ssGSEA scores for immune cells confirmed that the high-risk group had significantly fewer B cells, T cells, dendritic cells, macrophages, and nature killer cells than the low-risk group did (Figure 7A). ssGSEA score for immune functions demonstrated that pathways related to immune checkpoints, antigen-presenting cell, cell chemokine receptors, inflammation regulation, human leukocyte antigen, major histocompatibility complex, T cell functions, and type I/II interferon response were significantly downregulated in the high-risk group (Figure 7B).

The differences in the expression of 44 ICI-related genes between the high- and low-risk groups are shown in Figure 8A. The high-risk group had significantly lower expression levels of PD-1 and CTLA-4. Among patients whose disease expressed either of these ICI-related genes, low-risk patients were likely to have a better response to ICIs than high-risk patients were, through the Cancer Immune Atlas (TCIA) database (Figures 8B–E). These results suggest that a higher level of ferroptosis facilitated a higher immune response in the TME, leading to increase the efficacy of ICIs in CM patients. Then, we also used TIDE (Tumor Immune Dysfunction and Exclusion) score to evaluate the response to immunotherapy, and the results showed TIDE and Dysfunction scores were higher in low-risk group than high-risk group (Figures 8F,G). The prognostic performance in immunotherapy cohorts were displayed through survival cures, indicating that cohorts with low risk combined with high ICI score have the better prognosis than other cohorts, especially with PD-1 expression positive (Figure 8H).

The analysis demonstrated that the TMB had significantly difference between low- and high- risk groups (Figure 9A), and TMB has a significant impact on patient prognosis between low- and high- risk groups (Figures 9B,C). Then, we then identified the top 20 genes with the highest mutation rates in the two risk subgroups. The mutation rates of TTN, MUC16, BRAF, and DNAH5 were higher than 30% in both groups (Figures 9D,E). Mutation of the TTN gene was more common in the high-risk subgroup. Genetic mutations can affect the tumor response to chemotherapy and targeted therapy (Figure 9F); therefore, we investigated the association between the risk model and the efficacy of chemotherapy and targeted therapy drugs in patients with CM. We listed 10 common target drugs (Figure 9G) and 10 Chemotherapy drugs (Figure 9H) used for CM, such as sorafenib, Imatinib, ABT.888 (Veliparib), Docetaxel, and Paclitaxel (Figure 9H). Significant differences in the estimated IC50 between the two risk groups were observed, which suggest that the risk model might be used to identify potential biomarkers for chemotherapy and targeted therapy sensitivity.