The RNA-seq profiles of 471 CM and one non-tumor samples of the cancer genome atlas (TCGA),as well as 555 non-tumor skin specimens of the Genotype-Tissue Expression Project (GTEx), were obtained from the University of California Santa Cruz (UCSC). Batch normalization was performed in the data set by the sva R package. The differentially expressed genes (DEGs) were analyzed using Counts format profiles and converted into TPM format profiles for further analyses using the limma R package (Stupnikov et al., 2021). Their clinical data, such as overall survival (OS), and copy number, were also downloaded from UCSC. Furthermore, RNA-seq profiles, survival status, and OS time of CM patients were obtained from the GSE65904 and GSE54467 datasets for external validation from Gene Expression Omnibus (GEO). To reduce statistical bias, CM samples with missing OS values or short OS values (< 30 days) were excluded from all cohorts. As a result, 447 patients in the TCGA cohort and 278 in the external validation cohort were used for analyses. In addition, the fatty acid (FA) metabolism-related gene sets, including M14568, M22474, M34207, M34208, M29237, M23047, M15179, M16838, M25445, M23048, M13605, M11966, M15385, M16969, M14401, M34091, M13290, M11936, M14177, M12558, M16551, M13480, M15568, M18199, M26370, M26153, M40674, M18978, M26866, M26251, M5935, M36310, M37191, M11673, M29570, M27727, M27719, M14690, M22174, M6999, M13591, M10250, M15813, M25014, M16181, M12334, M17829, M25044, M15938, M23782, M22457, M40405, M40498, M699, M6995, M9927, M15531, M27854, M39440, and M39596, containing 745 fatty acid metabolism-related genes were extracted. Their intersection genes were identified in the TCGA, GTEx, and external validation cohorts using the VennDiagram R package. Moreover, the tumor mutation burden data of these cohorts were obtained from the TCGA database.

FA metabolism-related DEGs between the CM and non-tumor samples were identified, based on the criteria of a false discovery rate (FDR) < 0.05 and an absolute value Log2 fold change >1, using the limma R package in the TCGA and GTEx synthetic Counts matrix (Deng et al., 2021). This synthetic matrix had been normalized for identifying DEGs. Subsequently, 11 FA metabolism-related DEGs were identified using univariate Cox proportional hazard regression, and they had significant prognostic values (all p < 0.05) using the limma and survival R packages in the TCGA cohort (Zheng et al., 2021). These 11 FA metabolism-related DEGs with potent prognostic values were subjected to non-negative matrix factorization (NMF) analyses for sample clustering by NMF and survival R packages (Zhuo et al., 2020). The Kaplan–Meier survival curves of OS, t-distributed stochastic neighbor embedding (t-SNE), and Principal Component Analyses (PCA) of the clusters were analyzed using the survminer, Rtsne, ggplot2, and scatterplot3d R packages. Similar analyses were performed in the external validation cohort.

The copy number variation (CNV) frequencies of these FA metabolism-related prognostic genes were calculated, and their locations in human chromosomes were identified using the RCircos R package. For each prognostic gene, we compared the immune infiltration levels of CM patients with different somatic copy number alterations, such as deep deletion, arm-level deletion, and diploid/normal, using the Tumor Immune Estimation Resource (TIMER) (Kang et al., 2020).

The potential functions and pathways of each set of genes in the cluster were analyzed using the Gene Set Enrichment Analyses (GSEA) software (https://www.gsea-msigdb.org/gsea/login jsp) and Curated gene set (kegg. v7. 4. symbols. gmt). Multi-GSEA diagrams were drawn using the plyr, grid, and gridExtra R packages.

The immune infiltration of individual patients was analyzed using the CIBERSORT (R scrip v 1.03), and their immune scores, stromal scores, ESTIMATE (microenvironment) scores, and tumor purity were compared using the estimate R package, followed by visualizing them as the pheatmap using the ggtext packages (Gui et al., 2021). Subsequently, the immune cell infiltration, immune functions, and expression of genes for checkpoints were compared using the GSVA, GSEABase, ggpubr, reshape2, and limma R packages.

The tumor mutation burden (TMB) data in the “Masked Somatic Mutation” type were processed by VarScan2 and analyzed in clusters using the maftools package (Koboldt et al., 2012). Furthermore, their overall survival was estimated by K–M survival analysis, and a box plot was made using the survminer, survival, ggplot2, ggpubr, and ggExtra R packages.

According to the results of comprehensive immunogenomic analyses in The Cancer Immunome Atlas (TCIA), the response of each patient to immunotherapy was predicted and compared using the ggpubr R package (Charoentong et al., 2017). Additionally, the potential immunotherapeutic function of some drugs and their half-maximal inhibitory concentration (IC₅₀) in CM patients were predicted based on the data from Genomics of Drug Sensitivity in Cancer (GDSC) using the pRRophetic R package (Gui et al., 2021).

CM patients were randomly divided into training and testing sets. A risk model was generated using the least absolute shrinkage and selection operator (LASSO) in the glmnet R package. The risk score formula (Zheng et al., 2021): R i s k s c o r e = ∑ k = 1 n c o e f   ( g e n e   k )   *   e x p r   ( g e n e   k )

The coef (gene) meant the coefficient of the gene in the risk model, and expr (gene) was the expression of the gene in the risk model. The status, survival time, heatmaps, Kaplan–Meier survival analyses, and receiver operating characteristic (ROC) curves of individual patients in the training, testing, and entire sets were also analyzed using the pheatmap, survival, and timeROC R packages. For validation, a risk score of each patient in the external validation cohort was calculated, and their survival was estimated by the Kaplan–Meier survival analysis. The specificity and sensitivity of this risk model were evaluated by ROC analysis.

According to patient's clinical data and risk scores, the independent prognostic factors for worse survival were identified using univariate Cox and multivariate Cox regression analyses and used for the generation of a nomogram using the nomogramEx and nomogramFormula R packages. The concordance of the nomogram was analyzed by the ROC curves using the timeROC and rms R packages.

For predicting potential small-molecule drugs that might reverse high risk in CM, the whole overlapping genes, including the upregulated and downregulated genes, were submitted into the CMap database (https://portals.broadinstitute.org/cmap/). The drugs with enrichment scores between −1 and 0 were considered candidate drugs for CM (all p < 0.05) (Gui et al., 2021).

The tumor tissues were fixed in 10% formalin overnight and paraffin-embedded. The tissue sections (5 µm) were regularly stained with immunohistochemistry (IHC) staining. The stained tissue sections were photo-imaged and observed under a light microscope. The primary antibodies were used at 1:200 for CIDEA (Abcam, ab191193), ACOXL (Proteintech, 23366-1-AP), PLA2G2D (Abcam, ab47118), and 1:400 for KMO (Abcam, ab233529).

We made analyses with R version 3.6.3 (http://www.R-project.org) and its appropriate packages. All involved packages were described in MATERIALS AND METHODS. Data were analyzed with standard statistical tests as appropriate, while multiple testing was adjusted by the FDR method by R (Gui et al., 2021). Statistical significance was observed when p < 0.05.

A total of 745 FA metabolism-related genes were identified in CM from the GSEA (Supplementary Table S1). Similarly, 533 of them were identified in CM from the GTEx, TCGA, GSE65904, and GSE54467 (Figure 1A). There were 51 DEGs in the TCGA cohort. Of them, 17 were upregulated and the others were downregulated (Figure 1B). Further analyses indicated that 11 DEGs were associated with the prognosis of CM (all p < 0.05, Figure 1C). Of them, ALOX12B, CYP4F3, ALDH3A1, CIDEA, and BBOX1 were upregulated in CM, while the others were downregulated.

Based on comprehensive correlation coefficients and all heatmaps, we found the optimal total cluster number was set to k = 2 (Figures 1D,E). The heatmap (k = 2) indicated a clear boundary (Figure 1F). Compared with cluster 2, patients in cluster 1 had a worse OS (Figure 1G). To further verify the cluster distribution, the t-SNE, 2D PCA, and 3D PCA analyses clearly separated these prognostic genes (Figures 1H–J).

In the external verification cohort, the optimal total cluster number was k = 2 in NMF. The difference between the two clusters and their NMF ranks are displayed (Figure 2A, Supplementary Figures S1A,B). In addition, the t-SNE, PCA, and Kaplan–Meier survival curves of OS analyses revealed their distribution, and patients in cluster B had a worse OS in this population (Figures 2B–E).

Next, we calculated the CNV frequency of 11 prognostic genes and located them in human chromosomes. The percentages of gain CNV of KMO, KCNJ10, CYP4F3, and BBOX1 were higher than those of the loss, while the frequency of gain CNV of PLA2G2D, GPR31, IL4I1, ACOXL, ALDH3A1, ALOX12B, and CIDEA was lower than that of loss (Figure 2F). As a result, the KMO, KCNJ10, CYP4F3, and BBOX1 genes were marked in red, and the PLA2G2D, GPR31, IL4I1, ACOXL, ALDH3A1, ALOX12B, and CIDEA were marked in blue on RCircos 2D track plots (Figure 2G). Furthermore, we explored the changes in immune cell infiltration with a copy number alteration of prognostic genes in CM. Alterations of these genes were associated with infiltration levels of CD8⁺ T cells or other immune cells (Figure 2H). These indicated that clusters of CM with varying prognostic genes had different immune microenvironments, leading to different responses to immunotherapy (O'Sullivan et al., 2019).

We performed GSEA to explore the biological functions of these clusters. Cluster 2 and cluster 1 of genes were involved in the top 15 pathways (all p < 0.05, FDR <0.05, |NES|>1.8, Supplementary Figures S1C,D). In cluster 2, almost all of the enriched pathways were associated with immunity, such as chemokine signaling, and natural killer cell-mediated cytotoxicity (Figure 3A) (Chow and Luster, 2014; Taniguchi and Karin, 2018). On the contrary, six out of the top 15 pathways enriched by the genes in cluster 1 were associated with FA metabolism, such as arachidonic acid metabolism, four pathways were related to tumor growth, and two were related to drug metabolism (Figure 3B). According to the heatmap of immune cell infiltration, we felt that the cluster 2 had a higher immune cell infiltration level, lower tumor purity, and more active immunity (Figure 3C). In the immune cell bubble chart, cluster 2 of CM had more immune cell infiltrates (Figure 3D, Supplementary Data Sheet S1). Based on the single-sample GSEA scores for immune cells and immune functions, 13 types of immune cells, such as CD8⁺ T cells, and 12 immune functions, such as inflammation-promoting and check-point, had a higher score in cluster 2 of CM (ns means no significant difference, *p < 0.05, **p < 0.01, and ***p < 0.001, Figures 3E,F). Of course, almost all immune checkpoint genes, such as CD274 (PD-L1), CTLA4, HAVCR2 (TIME3),and LAG3, displayed a higher activation in cluster 2 (Figure 3G) (Galon and Bruni, 2019). Cluster 2 had a high immune score, stromal score, and ESTIMATE score (microenvironment score, Figure 3H). Therefore, cluster 2 of CM was considered the hot tumor to respond to immunotherapy because cluster 2 of CM displayed the characteristics of hot immune tumors, such as a high degree of CD8⁺ T cells, high immune score, more active immune functions, and higher expression of CTLA4, TIM3, and LAG3 (Galon and Bruni, 2019; Zheng et al., 2021). Previous studies have shown that a high TMB is significantly associated with improved OS and benefits from immunotherapy, such as CTLA-4 blockade, and the TMB has been considered a potential immunotherapy parameter (Snyder et al., 2014; Mahoney et al., 2015; Chan et al., 2019). Hence, we analyzed the TMB of CM in the TCGA cohort to explore the responses of these clusters of CM to immunotherapy. The TMB frequency in cluster 2 of CM (92.74%) was higher than that of cluster 1 (85.88%) (Figures 3I,J). Consistently, the high-TMB group (H-TMB) of CM patients had better OS than those in the low-TMB group (L-TMB) (Figure 3K). Stratification analysis revealed that CM patients with high TMB displayed a better OS than those with low-TMB in both clusters (Figure 3L).

Compared with cluster 1, CM in cluster 2 possessed significantly higher TMB (Figure 4A). Accordingly, we investigated immune checkpoint-related scores in these clusters of CM. Compared with cluster 1 CM, cluster 2 CMs had significantly higher both PD-1 and CTLA4 immunophenoscore (IPS, Figure 4B). In addition, CM patients in cluster 2 were likely to be more sensitive to 12 potential immunotherapy-related drugs with lower IC₅₀ than those in cluster 1 (Figure 4C) (Haikala et al., 2019; Song et al., 2019).

With the LASSO regression analyses, a risk model of seven genes was established after control of the first-rank value of Log(λ) at the minimum likelihood of deviance (Figures 5A,B). The risk score formula was Risk score = KMO × (−0.0032) + CYP4F3 × (−0.0054) + ACOXL × (−0.0220) + CIDEA × 0.0135 + IL4I1 × (−0.0039) + PLA2G2D × (−0.0033) + BBOX1 × 0.0048. In addition, analyses of the survival status, survival time, expression of seven genes, and OS of patients clearly separated between low- and high-risk groups of CM patients in the training, testing, and entire sets of cases (Figures 5C–N). Hence, CM patients in the low-risk group displayed a better prognosis. The AUC for 1-, 2- and 3-year OS in the training set of cases was 0.671, 0.703, and 0.701, that of the testing set of cases was 0.664, 0.665, and 0.654, and that of the entire sets of cases was 0.663, 0.682, 0.676, respectively (Figures 5O–Q). After analyzing the clinical characteristic of the patients, we found our risk model was also available in age, gender, tumor stage, T stage. and N stage (Supplementary Figure S2A). Furthermore, stratification of CM patients was carried out into low- and high-risk groups, and the t-SNE, 2D PCA, and 3D PCA clearly separated them (Supplementary Figures S3A). Further validation revealed that CM patients in the low-risk group exhibited significantly better OS than those in the high-risk group (p = 0.002, Figure 5R). The AUC for 1-, 2-, and 3-year OS in the external validation cohort were 0.668, 0.684, and 0.690, respectively (Figure 5S).

The univariate Cox (uni-Cox) and multivariate Cox (multi-Cox) regression analyses indicated that patient’s age, tumor T stage, tumor N stage, and risk score were independent risk factors for prognosis of CM patients. Their independent risk factors were associated with a worse prognosis for CM patients with hazard ratios and a 95% confidence interval (CI) (Figures 6A,B). With these four independent prognostic indexes, a nomogram was generated (Figure 6C). The AUC of a nomogram for 1-, 2-, and 3-year OS of CM patients was 0.766, 0.816, and 0.816, respectively (Figure 6D). The calibration plots of the nomogram for 1-, 2-, and 3-year OS of CM patients had a good concordance with the prediction (Figure 6E).

Among the top 15 signaling pathways, 11 signaling pathways were associated with immunity, such as the T cell receptor signaling (all p < 0.05, FDR <0.05, |NES|>2.3) (Supplementary Figure S3B) (Courtney et al., 2018). The immune heatmap exhibited that CM in the low-risk group had more types of immune cells (Supplementary Figure S3C, Supplementary Data Sheet S2). The levels of CD8⁺ T, B cell infiltrates, and other types of immune cells were also correlated with lower risk scores (Supplementary Figures S3D). Furthermore, CM in the low-risk group had higher immune scores (Supplementary Figures S3E). Thus, we tried to explore potential immunotherapy-related drugs for these risk groups of CM. As a result, we found that four drugs, such as metformin, had significantly different IC₅₀ between both risk groups of CM (Figure 7A) (O'Sullivan et al., 2019; Cha et al., 2018). A total of five most related small-molecule drugs, ambroxol, tiletamine, mimosine, esculetin, and pizotifen, were identified (all p < 0.05), based on upregulated and downregulated gene expression between the low- and high-risk groups of CM (Supplementary Data Sheet S3) (Gui et al., 2021). Their 2D and 3D structure tomography are is in Figures 7B,C.

Finally, we further validated the expression of these key genes in normal skin and CM tissues. Melanoma arises from melanocytes in the epidermal layer of the skin, so we focused on the expression of these proteins in the epidermis. First, we searched the immunohistochemical slide information through the public database on the HPA website and found that the expression of IL4I1 in normal skin epidermis was higher than that in CM tissue, while BBOX1 protein was not detected in normal skin (Figures 8A,B). Next, our own immunohistochemical results showed that the expression of PLA2G2D, ACOXL, and KMO was higher in normal epidermal tissue than in CM tissue; however, the expression of CIDEA was reversed (Figures 8C–F).