Transcriptome data related to SKCM were retrieved from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). A total of 477 SKCM samples were obtained from the TCGA dataset. After samples lacking survival information were removed, 457 samples with complete survival information were retained for subsequent analysis. Clinical information samples, including age, gender, tumor stage, M, N, and T staging, were extracted. After samples lacking clinical information were deleted, 236 samples remained for subsequent clinical information related analyses.

Additionally, datasets GSE46517 and GSE65904 were obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). The GSE46517 dataset contained 121 samples, which included metastatic melanoma, primary melanoma, nevus, and normal skin samples. The metastatic melanoma, melanoma, and nevus samples were excluded, and 31 primary SKCM samples along with 7 normal samples were retained for analysis. The GSE65904 dataset contained 210 SKCM samples, and all of them were included in the subsequent analysis. A total of 19 PANoptosis-related genes (PRGs) were sourced from existing literature ( Supplementary Table 1 ) (15).

In TCGA-SKCM dataset, the enrichment scores of the PRGs in the SKCM samples were calculated using the ssGSEA method in the GSVA package (version 1.42.0) of R language. The parameters were set as gsva(as.matrix(exp), geneSet, method=‘ssgsea’, kcdf=‘Gaussian’, abs.ranking=TRUE). Subsequently, the enrichment scores of PRGs were defined as the PRGs score. To explore the impact of the PANoptosis score on the survival of SKCM patients, SKCM patients were divided into a high-PRG-score cohort and a low-PRG-score cohort based on the median PANoptosis score (16). This approach effectively balanced the sample sizes and provided a stable statistical foundation for subsequent survival analysis. Kaplan-Meier (K-M) survival analysis was then performed to compare survival rates between the two groups (p < 0.05) (17). Differentially expressed genes (DEGs) between the two groups were identified using the DEseq2 package in R (|log₂FC| > 1 and P adj. < 0.05) (18) Volcano and heat maps were generated using the ggplot2 and pheatmap packages in R, respectively, to visualize the results (19).

To identify key PANoptosis-associated module genes in SKCM samples, weighted gene co-expression network analysis (WGCNA) was conducted. Firstly, outlier samples in the TCGA-SKCM dataset were removed through cluster analysis. Then, the coefficient of determination (R^{^2}) of the scale-free network was set to 0.85, and the optimal soft threshold was obtained when (R^{^2}) first exceeded the critical value of 0.85 and the average connectivity of the co-expression network was close to 0 (20). Next, according to the screened soft threshold, a scale-free network was constructed, and the minimum number of genes in each module was set to 200, so that genes could be effectively grouped into multiple modules. Subsequently, the PANoptosis score was taken as the trait. Through Pearson correlation analysis of the correlation between the module and the trait, the module with the strongest significant correlation was selected as the key module.²⁰ Genes with |Module Membership (MM)| > 0.8 and |Gene Significance (GS)| > 0.5 in the key module were considered as key module genes (21, 22).

To identify DEGs in samples from normal and SKCM cohorts, the limma R package was used with the GSE46517 dataset (|log₂FC| > 1 and P adj. < 0.05) (23). A volcano plot was created using the ggplot2 R package (version 3.3.6) to illustrate the filtering results (19). Heatmaps of DEG expression were visualized using the pheatmap R package (version 1.0.12) (23). Differentially expressed PANoptosis-related genes (DEPRGs) were identified by intersecting DEGs between normal and SKCM samples, key module genes, and DEGs from high and low-PANoptosis score groups.

Finally, Kyoto Encyclopedia of Genes and Genomes and Gene Ontology (GO) analyses for the DEPRGs were performed using the clusterProfiler R package (version 4.4.4) (24).

In the TCGA-SKCM dataset, 457 samples were divided into a validation set (137 samples) and a training dataset (320 samples). Univariate Cox analysis was performed on the DEPRGs in the training dataset to identify prognostic-related genes (HR ≠ 1, p < 0.05). Subsequently, the least absolute shrinkage and selection operator (Lasso) was applied using the R package glmnet. Based on 10-fold cross-validation, prognostic-related genes with the smallest cross-validation error value (λ value) and non-zero regression coefficients were selected as the characteristic genes for SKCM patients. A multifactorial Cox analysis was then conducted on these characteristic genes to identify prognostic genes and build a prognostic model for SKCM (25). The variance inflation factor (VIF) values of the multivariate Cox regression model were calculated using the vif function, and the model performance metric, the C-index, was computed (C-index > 0.6). Subsequently, the proportional hazards (PH) assumption of the multivariate Cox regression model was tested using cox.zph (p > 0.05). To investigate whether severe multicollinearity existed among the prognostic genes, a linear model was constructed using the `lm` function with each prognostic gene as the dependent variable and the other genes as independent variables. The VIF for each gene was then calculated using the vif() function from the `car` package. A VIF value less than 5 indicated the absence of severe multicollinearity.

To evaluate the prognostic value of the model, patients with SKCM in the training dataset were classified into high- and low-risk groups based on the median risk score (25). K-M survival analysis was performed to compare survival rates between the two groups. Receiver operating characteristic (ROC) curves for the prognostic model were generated using the survival ROC R package (17). The robustness of the model was further validated using the external GSE65904 dataset and the TCGA-SKCM validation set.

The independent prognostic value of clinicopathological parameters (including age, history of neoadjuvant treatment, prior treatment.diagnoses, prior systemic therapy, M stage, N stage, T stage, gender and tumor_stage) and risk scores was assessed using univariate and multivariate Cox proportional hazards regression analyses in the training set (26, 27). A nomogram model was then constructed by integrating factors with independent prognostic significance (17). The rms package (version 6.5.1) was used to generate a calibration curve, which was applied to assess the precision of the nomogram’s prognostic ability (28).

To explore the biological functions and pathways associated with the prognostic genes in SKCM, GSEA was performed using TCGA-SKCM data. The DESeq2 R package was used to identify differential gene expression between the high- and low-risk groups in the training dataset. These genes were ranked according to their log₂ fold-change (log₂FC) values and subjected to GSEA through the clusterProfiler and org.Hs.eg.db R packages (24).

To examine differences in immune cell infiltration between the high- and low-risk cohorts, the CIBERSORT algorithm was applied to determine the distribution of 22 immune cell types in TCGA-SKCM (p < 0.05). A bar plot, created with the ggplot2 R package, was used to visualize the immune cell expression in the two risk groups (25). The Wilcoxon test was employed to compare the distribution differences of immune cells between the groups (p < 0.05) (29). Finally, the psych R package (version 2.3.12) was used to explore the relationships between differentially infiltrated immune cells and prognostic genes through Spearman correlation analysis (30).

The molecular regulatory mechanisms of prognostic genes were explored in future studies. The miRNAs associated with prognostic genes were identified using the miRTarBase database, and the upstream miRNAs with a connectivity greater than 1 and supported by experimental results were screened (31). Next, miRNA-associated long non-coding RNAs (lncRNAs) were predicted through the starBase database, applying a screening criterion of clipExpNum > 10 (27). Competitive endogenous RNA (ceRNA) networks were constructed using Cytoscape (31).

To explore drugs associated with prognostic genes and identify potential treatments for SKCM, potential therapeutic agents for SKCM were predicted using the Drug-Gene Interaction Database (DGIdb) (http://dgidb.genome.wustl.edu/), PubChem (https://pubchem.ncbi.nlm.nih.gov/), and the Therapeutic Target Database (TTD) (https://db.idrblab.net/ttd/), based on the prognostic genes (32–34). The drug-prognostic gene network was then visualized using Cytoscape.

In this study, RT-qPCR analysis was performed on HSF, MV3, SK-MEL-28, and WM-115 cell lines purchased from the ATCC database when cell confluence reached 70%-90%. Three samples were selected for each cell type. Total RNA was extracted from these cells using TRIzol reagent (Ambion, Austin, USA), and incubated for 15 minutes. One microliter of RNA solution was measured using the NanoPhotometer N50. Reverse transcription to synthesize complementary DNA (cDNA) was performed using the SureScript First-strand cDNA Synthesis Kit (Servicebio, Wuhan, China). Following cDNA synthesis, the solution was diluted 5 to 20 times with RNase/DNase-free distilled water (ddH₂O). A 10 μl reaction system was configured for qPCR, which contained 3 μl of cDNA, 5 μl of 2xUniversal Blue SYBR Green qPCR Master Mix, 1 μl of 10 μM Forward primer, and 1 μl of 10 μM Reverse primer. RT-qPCR amplification was carried out using a CFX96 Real-Time Fluorescence Quantitative PCR Instrument (Bio-Rad, California, USA) for a total of 40 cycles. Finally, Ct values of prognostic genes and calibrator gene Glyceraldehyde-3-Phosphate Dehydrogenase(GAPDH) were employed to compute relative expression of every prognostic gene in each sample by 2^−△△Ct method. The t-test of Graphpad Prism 5 (v 8.0) software (35). GAPDH was used to assess RT-qPCR data (p < 0.05). Detailed information about the RT-qPCR reaction system, primer sequences, and amplification conditions is provided in Supplementary Tables 2–4 .

Statistical analyses were performed using R software (v 4.1.0), and the Wilcoxon test was used to compare pairs within cohorts. A p-value of < 0.05 was considered statistically significant.

This study, titled The role and prognostic value of PANoptosis-related genes in skin cutaneous melanoma, has been approved by the ethics committee of Shanxi Bethune Hospital. The approval number and date of approval are as follows: [YXLL-2023-087] and [2023-4-11].

In the TCGA-SKCM dataset, individuals in the low PANoptosis score group exhibited significantly lower survival probabilities compared to the high PANoptosis score group ( Figure 1A ). Differential expression analysis revealed 1460 DEGs between the high and low PANoptosis score groups (|log₂FC| > 1 and P adj. < 0.05), including 219 down-regulated genes and 1241 up-regulated genes ( Figures 1B, C ). All samples in the TCGA-SKCM dataset were successfully clustered, with no outliers observed ( Figure 1D ). To construct a co-expression network and identify modules, the soft threshold power of seven was selected to compute adjacencies ( Figure 1E ). Dynamic tree-cutting identified eight modules ( Figure 1F ), with the MEblue module showing the strongest correlation with PANoptosis scores (cor = 0.83, p < 0.001) ( Figure 1G ). Of the 3801 genes in the MEblue module, 426 were selected as key module genes ( Figure 1H ).

Variable analysis of the GSE46517 dataset identified 854 DEGs between normal and SKCM samples (P adj. < 0.05 and |log₂FC| > 1), comprising 472 down-regulated genes and 382 up-regulated genes ( Figures 2A, B ).

A total of 26 DEPRGs were identified by intersecting 854 DEGs from normal vs. SKCM samples, 1460 DEGs from high and low PANoptosis score groups, and 426 key module genes ( Figure 2C ). GO analysis revealed that DEPRGs were significantly enriched in lymphocyte differentiation, monocyte differentiation and migration, cytokine activity, and chemokine activity ( Figure 2D ).

KEGG analysis revealed significant enrichment of DEPRGs in cytokine and chemokine signaling pathways, as well as their corresponding receptors ( Figure 2E ).

Univariate Cox analysis identified 26 prognosis-related genes in the training dataset (HR ≠ 1, p-value < 0.05) ( Table 1 ). Following this, LASSO analysis showed that when the minimum λ value (lambda.min = 0.011) was reached, the minimum error was obtained, and 11 characteristic genes with non-zero regression coefficients were derived ( Figures 3A, B , Table 2 ). Multifactorial Cox analysis further identified seven prognostic genes and constructed a prognostic model ( Figure 3C ). The VIF value of each gene was less than 5 ( Table 3 ), and the C-index of the model was 0.64. The PH assumption test showed that the p-values of all genes were greater than 0.05 ( Table 4 ), indicating that these genes had a certain predictive ability. The seven prognostic genes were CD8A, ADAMDEC1, CD69, CRIP1, LSP1, BCL11B, and CCR7. The multicollinearity test for the prognostic genes showed that the VIF value for each gene was less than 5, indicating that there was no significant multicollinearity ( Supplementary Figure 1 ). Patients in the high-risk cohort, who exhibited shorter survival times, were predominantly found in the DEPRGs cases ( Figure 3D ). The high-risk group demonstrated significantly worse survival outcomes (P < 0.0001) ( Figure 3E ). The prognostic model proved to be reliable in assessing DEPRGs in the training dataset ( Figure 3F ).

Similar results were obtained in both the GSE65904 dataset and the validation set, with patients in the high-risk group showing slower survival ( Figures 4A, B ). The 1-year, 3-year, and 5-year ROC curves for SKCM patients were plotted. The AUC values were found to be 0.646, 0.657, and 0.681 in GSE65904 set, respectively, and 0.625, 0.660, and 0.646 in validation set, respectively ( Figures 4C, D ). The above results demonstrate that the prognostic model is reliable in evaluating SKCM.

In the training set, univariate Cox analysis revealed significant associations between risk score, age, TNM stage, and clinical stage with the survival of patients with SKCM (P < 0.05). While, history of neoadjuvant treatment, prior treatment diagnoses, prior systemic therapy, and gender were found to have no significant impact on the survival of SKCM (P > 0.05) ( Figure 5A ). Multivariate Cox analysis identified N stage, T stage, and risk score as independent prognostic factors (P < 0.05) ( Figure 5B ).

A nomogram incorporating T stage, N stage, and risk score was constructed to reliably predict the prognosis of patients with SKCM ( Figures 5C, D ).

KEGG analysis results indicated that genes in both groups were predominantly associated with interactions between cytokines and cytokine receptors, intestinal immunity networks for IGA production, systemic lupus erythematosus, and allogeneic transplant rejection, among others ( Figure 6A ). GO biological process analysis revealed substantial enrichment of genes in both groups in processes such as adaptive immune response, epidermal and keratinocyte differentiation, and biological processes involving intermediate filaments ( Figure 6B ). GO cellular component analysis showed significant enrichment in components like the cornified envelope, immunoglobulin complex, and intermediate filaments ( Figure 6C ). GO molecular function analysis identified key functions such as antigen binding, immune receptor activity, immunoglobulin receptor binding, and MHC protein complex binding ( Figure 6D ). According to the oncogenic signature database, GSEA highlighted a primary association with the upregulation of RPS14, SNF5, and STK33 ( Figure 6E ).The above results established associations between SKCM prognostic genes and key biological processes such as immune regulation, cell differentiation, and oncogenic pathways through systematic functional genomics analysis. This not only provided a logical framework for explaining the molecular mechanisms of prognostic differences but also laid a data foundation for developing combined targets for immunotherapy.

The distribution of immune cells was analyzed in samples from the high- and low-risk groups ( Figure 7A ). Significant differences were observed in the distribution of nine immune cell types between the two groups ( Figure 7B ). Prognostic genes showed an inverse correlation with Macrophage M0 and a positive correlation with most differentially infiltrating immune cells ( Figure 7C ). Through immune cell infiltration analysis, the dynamic associations between SKCM prognostic genes and the tumor immune microenvironment were revealed, providing key insights for deciphering the immunological mechanisms of prognostic differences, developing immunotherapy biomarkers, and optimizing combination treatment strategies.

miRNA prediction was performed for three prognostic genes (CD69, LSP1, and BCL11B), leading to the construction of a ceRNA network comprising three mRNAs, 56 miRNAs, and 27 lncRNAs ( Figure 8A ). Among these, SNHG3, AL355075.4, and LRRC75A-AS1 regulate BCL11B expression by modulating hsa-miR-326 and hsa-miR-330-5p. The non-coding RNA regulatory mechanisms of SKCM prognostic genes were revealed above, and potential ceRNA regulatory axes and therapeutic targets were screened out, providing a theoretical basis for deepening the understanding of tumor molecular mechanisms and developing novel biomarkers or RNA therapies.

Seven potential therapeutic agents for SKCM were identified based on CD8A, LSP1, CD69, and CCR7, while no agents were predicted based on ADAMDEC1, CRIP1, or BCL11B ( Figure 8B ). The identified therapeutic agents included Fluoropeptide vaccine, LIPO-4, Anti-VEGFR2 CD8 cell therapy, MPPG, DT-1687, TOCILIZUMAB, and JBH492. By integrating prognostic gene expression and drug sensitivity data, the screening range of therapeutic agents for SKCM was narrowed, providing candidate regimens for personalized treatment.

As presented in Figure 9 , RT-qPCR analysis demonstrated that in SKCM cell lines (MV3, SK-MEL-28, and WM-115), expression levels of CD8A, ADAMDEC1, CD69, and CRIP1 were significantly lower compared to the HSF cell line, while BCL11B expression was significantly higher. No significant difference in the expression of LSP1 and CCR7 was observed.