In this study, a total of four bulk sequencing datasets for cutaneous melanoma were downloaded and utilized. These datasets include The Cancer Genome Atlas (TCGA) cohort, GSE19234 cohort, GSE22153 cohort, and GSE65904 cohort (Bogunovic et al., 2009; Jönsson et al., 2010; Cabrita et al., 2020). The “sva” method was applied to standardize the downloaded bulk data, effectively mitigating any noticeable batch effects. Additionally, we matched gene expression data with clinical survival data. Log2 transformation was uniformly applied to all the data.

The single-cell sequencing data for cutaneous melanoma (GSE115978) was obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) in count format (Jerby-Arnon et al., 2018). To standardize this data, both the “LogNormalize” and “vst” methods were applied. High variable genes were selected using the “FindVariableFeatures” function. Subsequently, PCA (Principal Component Analysis) was performed for dimension reduction, and the batch effects were mitigated using the harmony function. Cell types were annotated based on highly expressed markers within each cluster, and visualization was carried out using t-SNE and UMAP techniques.

Inflammation-related genes were sourced from the Genecards database (https://www.genecards.org). Using the search bar, “inflammation” was entered as the keyword for retrieval. The top 100 genes ranked by their scores were exported (Supplementary Table S1). These genes are considered highly relevant to inflammation and will be used for subsequent analysis.

The gene expression score was initially computed for the melanoma single-cell sequencing data GSE115978 using the “AUCell” R package. Following this, the 100 genes related to inflammation were extracted to derive inflammation-related area under the curve (AUC) scores. These scores were then divided into two groups, namely, the high-AUC group and the low-AUC group, based on the median value. Differential gene expression analysis was subsequently conducted to identify genes that exhibited differential expression between these two groups. To explore the differential enrichment pathways between the two groups, Gene set variation analysis (GSVA) analysis was employed.

Once the individual cell inflammatory AUC scores were acquired, the “cor.test” function was employed to perform Spearman correlation analysis to identify genes correlated with these AUC scores. Subsequently, the first 150 genes found to be associated with the AUC scores were intersected with the differentially expressed genes identified between the two AUC groups. This final outcome was regarded as the most pertinent set of genes related to inflammation in the single-cell sequencing data.

In this study, we conducted a Weighted Gene Co-expression Network Analysis (WGCNA) of a set of 100 inflammation-related genes on TCGA melanoma dataset. Following data normalization, we constructed a co-expression network using the “WGCNA” R package, identified gene modules, and assessed their correlation with inflammation-related traits. Hub genes within modules were determined based on intramodular connectivity, and functional annotation was performed to elucidate the biological processes associated with these modules. This analysis aimed to uncover co-expression patterns, potential key regulators, and functional insights related to inflammation in the TCGA dataset. Finally, we intersected the key genes in the key modules obtained by WGCNA with the key genes in the single-cell sequencing data, and the remaining genes were used to construct a prognostic model.

In this study, we employed a multi-step approach to construct prognostic models using a selected set of genes derived from unifactor COX regression analysis on the TCGA dataset and applied machine learning techniques in three independent GEO cohorts. First, unifactor COX regression analysis was conducted on the TCGA dataset to identify genes strongly associated with prognosis. Subsequently, a variety of machine learning algorithms, including RSF, COXBoost, Enet, GBM, Lasso, plsRcox, Ridge, StepCox, and Survivor-SVM, were utilized individually and in combinations to construct prognostic models (Wang Q. et al., 2023; Wang D. et al., 2023; Pei et al., 2023). During the construction of the model, the TCGA cohort is used as the training set and the three GEO cohort is used as the validation sets. The Concordance Index (C-index) comparison of the models is performed in validation sets. To assess the performance of these models, the C-index was calculated for each model, ranking them from high to low. The “ComplexHeatmap” package in R was used to visualize the results, providing a comprehensive overview of the prognostic model performance across different algorithms and their combinations. This analysis aimed to identify the most effective models for predicting prognosis in the context of the three GEO cohorts, with the C-index serving as a key metric for evaluation.

The prognostic model was applied to multiple cohorts, including the TCGA cohort, GSE19234, GSE22153, and GSE65904. For each patient in these cohorts, a risk score was computed. Subsequently, patients were categorized into high-risk and low-risk groups based on the median value derived from these four cohorts. To assess survival outcomes, we utilized the Kaplan-Meier method to generate survival curves and conducted log-rank tests to establish statistical significance. Furthermore, we developed a nomogram that integrates the risk scores with various clinical characteristics to predict survival probabilities at 1, 3, and 5 years. To gauge the nomogram’s accuracy, we constructed calibration curves. Additionally, we evaluated the model’s robustness by generating ROC curves and calculating the area under the curve (AUC).

To gain deeper insights into the prognostic model’s relevance in the tumor microenvironment, we conducted immune and mutational analyses, along with predicting differences in immunotherapy response between the high-risk and low-risk groups. These analyses were carried out within the TCGA cohort. Initially, we divided patients in the TCGA cohort into high-risk and low-risk groups based on their risk scores, which were previously calculated. Next, we obtained data regarding immune cell infiltration and matched patient groups with their respective immune cell infiltration profiles. The results were then visualized using the “pheatmap” R package. Similarly, we calculated the tumor mutation load (TMB) in a manner analogous to the immune cell infiltration analysis. TMB data from the TCGA cohort were linked with the risk score groups to identify genes with high mutation rates. Visualization was accomplished using the “ComplexHeatmap” package. Furthermore, we utilized the “estimate” R package to compute immune scores, stromal scores, Estimate scores, and tumor purity for both the high-risk and low-risk groups. To evaluate the potential efficacy of immunotherapy, we employed immunotherapy data from the TCIA database to predict TCIA scores for the high-risk and low-risk groups, thereby assessing differences in immunotherapy response.

The human malignant melanoma cell lines, namely, A375 and WM115, were sourced from the cell repository established by the Type Culture Collection of the Chinese Academy of Sciences. These cell lines were cultured in DMEM (Gibco BRL, Rockville, MD, United States), complemented with 10% fetal bovine serum (FBS, Gibco BRL, Rockville, MD, United States), and a 1% penicillin-streptomycin solution. Incubation was carried out at 37°C in a 5% CO₂ environment.

To conduct gene knockdown experiments, small interfering RNA (siRNA) and its corresponding negative control (NC) were custom-synthesized on our behalf by Ribobio, based in Guangzhou, China. The specific siRNA sequences targeting STAT1 can be found in Supplementary Table S2. Subsequently, we performed transfections on A375 and WM115 cells using Lipofectamine 3000, following the manufacturer’s instructions provided by Invitrogen, situated in Carlsbad, CA, United States.

The evaluation of cell proliferation in both A375 and WM115 cell lines involved the utilization of Cell Counting Kit 8 (CCK-8). Following transfection, cells were seeded into 96-well plates overnight, with each well containing 2000 cells (maintained at 37°C under 5% CO₂). Subsequently, a mixture comprising 10 µL of CCK-8 labeled reagent (A311-01, Vazyme, Nanjing, China) and 90 µL of serum-free medium was added to the cells. They were then incubated in darkness at 37°C for 2 h, and the absorbance at 450 nm was measured using an enzyme marker (A33978, Thermo, Waltham, MA, United States).

For the clony formation assays, transfected A375 and WM115 cells were cultured in 6-well plates for approximately 12 days. Following this incubation period, cells were subjected to staining with 0.1% crystal violet for 30 min, followed by a PBS wash. The colonies exceeding 1 mm in diameter were then enumerated.

Invasion and migration assays were executed by employing a 24-well transwell chamber (BD Biosciences, Franklin Lakes, NJ, United States). In the upper chamber, A375 cells were suspended in 200 µL of serum-free medium, while the lower chamber contained 600 µL of medium supplemented with 10% fetal bovine serum. After incubating for 48 h, the cells that had penetrated the upper membrane surface were meticulously removed. The remaining cells were fixed with 4% paraformaldehyde and subjected to staining with 0.1% crystal violet (Solarbio, Beijing, China). Following this, three distinct fields of view were observed and quantified under a light microscope at ×200 magnification.

In R software (version 4.3.2), the analysis involved several steps. Single-cell analysis was executed using the Seurat package, while the assessment of gene heterogeneity associated with inflammation was conducted through the AUCell method. Highly correlated genes were identified using Spearman correlation analysis, and modules linked to inflammation were pinpointed using the WGCNA method. Prognostic models were constructed using machine learning techniques such as RSF, COXBoost, Enet, GBM, Lasso, plsRcox, Ridge, StepCox, and survival-SVM. Survival analysis was performed using the Kaplan-Meier (K-M) method, and the robustness of the models was assessed using ROC curves. Throughout the study, all experiments were rigorously repeated three times to ensure the reliability of the findings, and statistical significance was determined with a threshold of p < 0.05.

Following harmony batch removal and tSNE dimensionality reduction clustering, the entire cell population was segregated into 19 clusters, as depicted in Figure 1A. Subsequently, through the annotation of characteristic genes within these clusters, they were identified as B cells, cancer-associated fibroblasts (CAFs), endothelial cells, macrophages, malignant (Mal) cells, NK cells, and T cells, as illustrated in Figure 1B. To assess the activation of inflammation-related genes at the single-cell level, AUCell analysis was employed, classifying cells into either the high-AUC group or the low-AUC group based on the median value. Notably, NK cells, T cells, and macrophages exhibited elevated inflammation-related scores. Further analysis involved the extraction of differentially expressed genes between the high-AUC and low-AUC groups, as depicted in Figure 1C. Spearman correlation analysis was conducted to pinpoint the top 150 genes strongly associated with the inflammatory AUC score, as showcased in Figures 1D, E. Subsequently, 67 genes were selected for further analysis through the intersection with the pool of differentially expressed genes. The identification of differentially activated pathways between the high-AUC and low-AUC groups was accomplished using GSVA analysis, encompassing pathways such as allograft rejection and inflammatory response, as presented in Figure 1F.

WGCNA analysis of cutaneous melanoma revealed the clustering of various genes into distinct modules, and their relationships with inflammatory phenotypes were assessed. Notably, the blue module exhibited the most robust positive correlation with inflammatory phenotypes, as illustrated in Figures 2A–C. Consequently, a total of 3,162 genes within the blue module were identified by WGCNA as highly associated with inflammation. Subsequently, the 3,162 genes obtained from WGCNA were intersected with the 67 genes acquired through single-cell sequencing analysis, resulting in the identification of 46 genes for the construction of machine learning models (Figure 2D).

Following the acquisition of the aforementioned 46 genes, the identification of prognostic-related genes commenced with univariate COX regression analysis. Subsequently, a prognostic model was constructed using a combination of machine learning algorithms, including RSF, COXBoost, Enet, GBM, Lasso, plsRcox, Ridge, StepCox, and survival-SVM. After thorough evaluation, it was determined that the CoxBoost + GBM method yielded the most significant prognostic model, with a C index of 0.604 (Figure 3). This prognostic model consists of six genes, namely, FCGR2A, TLR2, STAT1, IFIT2, LAP3, and PARP14.

In the TCGA cohort, GSE19234 cohort, GSE22153 cohort, and GSE65904 cohort, the application of this prognostic model enabled the stratification of all patients into high-risk and low-risk groups. Notably, across all four cohorts, a distinctly poorer prognosis was consistently observed in the high-risk group, as depicted in Figures 4A–D (p < 0.05).

Utilizing TCGA data, a nomogram was developed to assess the 1, 3, and 5-year survival probabilities of patients by incorporating their risk scores, gender, age, and stage. The nomogram demonstrated that the estimated survival rates for patients at 1, 3, and 5 years were 0.987, 0.919, and 0.863, respectively, as illustrated in Figure 5A. The calibration curve further confirmed the accurate performance of the nomogram, as depicted in Figure 5B.

Subsequently, ROC curves were generated for 1, 3, 5, and 7-year survival predictions within the TCGA cohort. It was observed that both the risk score and the nomogram exhibited higher AUC values compared to other clinical indicators, affirming the robust predictive capabilities of the prognostic model and the nomogram, as depicted in Figure 6A. Furthermore, time-ROC curves spanning 1, 3, 5, 7, and 10 years were constructed to assess the accuracy of the prognostic model in the TCGA cohort, GSE19234 cohort, GSE22153 cohort, and GSE65904 cohort, demonstrating the model’s effectiveness, as illustrated in Figure 6B.

In our analysis, we initially constructed the immune cell infiltration landscape, distinguishing between the high-risk group (left) and the low-risk group (right) as shown in Figure 7A. Notably, it became evident that the low-risk group potentially exhibited a higher level of immune cell infiltration. Further exploration of the differences in Tumor Mutational Burden (TMB) between these two groups unveiled genes with the highest mutation frequency, as illustrated in Figure 7B. Investigating immune checkpoint-related genes, we observed that the low-risk group displayed elevated expression levels of immune checkpoint-related genes, suggesting a potential advantage in terms of immunotherapy responsiveness, as indicated in Figure 7C. Additionally, utilizing the TICA score, we found that the low-risk group had a higher TCIA score, potentially indicating a more favorable immunotherapy response, as depicted in Figure 7D. Various scores related to the tumor microenvironment revealed that the risk score exhibited a negative correlation with stromal score, immune score, and ESTIMATE score, while displaying a positive correlation with tumor purity (Figure 7E) (*p < 0.05, ***p < 0.001).

By constructing a Protein-Protein Interaction (PPI) network, our analysis revealed that among the six model genes, STAT1 occupied a central position within the network, suggesting it may have a crucial hub role, as depicted in Figure 8. Consequently, we chose to validate the function of STAT1 in cutaneous melanoma through cellular experiments. First, the Cell Counting Kit-8 (CCK-8) experiment demonstrated a significant decrease in the activity of A375 and WM-115 cell lines following STAT1 knockdown, as shown in Figure 9A(***p < 0.001). Furthermore, clonal formation assays indicated a marked reduction in colony formation ability for both A375 and WM-115 cell lines after STAT1 knockdown, as illustrated in Figures 9B, C (**p < 0.01, ***p < 0.001). Transwell experiments revealed a significant decrease in the migration ability of A375 and WM-115 cell lines upon STAT1 knockdown, as depicted in Figures 9D, E (***p < 0.001).