Multi-Omics Integration with Machine Learning and Molecular Docking Reveals Crosstalk Mechanisms and Drug Candidates in Metastatic Melanoma and Vitiligo

Introduction

Melanoma is a malignant tumor of the skin that arises from the malignant proliferation of melanocytes. It is characterized by strong aggressiveness, early metastasis and a poor prognosis.¹ The American Society of Clinical Oncology (ASCO) has observed a notable increase in the incidence and lethality of melanoma over recent decades, with the disease now representing a significant public health burden.² It is widely accepted that sun exposure and ultraviolet radiation (UVR) are significant contributing factors in the development of melanoma.³ The advent of immune checkpoint inhibitors in conjunction with conventional surgical treatments has enhanced the survival prospects of melanoma patients. However, a proportion of patients remain unresponsive or exhibit drug resistance.⁴

In clinical practice, an association has been identified between vitiligo and melanoma. A nonconcurrent cohort demonstrates a significantly reduced incidence of melanoma in patients with vitiligo, despite the potential for increased exposure to UVR.⁵ Additionally, another study finds that individuals with vitiligo exhibited a threefold reduction in the likelihood of developing melanoma.⁶ Meanwhile, a meta-analysis indicates that patients with stage III–IV melanoma who have developed vitiligo-like depigmentation as a result of immunotherapy exhibit a two- to four-fold reduction in the risk of disease progression and mortality, respectively, in comparison to those who have not developed this phenomenon.⁷ These findings suggest that there may be a protective effect associated with vitiligo in the development of melanoma. Vitiligo is an autoimmune disease characterized by the destruction of melanocytes, which results in the loss of pigmentation in the skin.⁸ This mechanism is diametrically opposed to that observed in melanoma, a malignant proliferation of melanocytes. It may represent a potential common mechanism of action for vitiligo and melanoma.

The current research indicates that the correlation between vitiligo and melanoma is believed to be the consequence of an immune response against melanoma and normal melanocytes.⁹ For example, autoantibodies isolated from patients with vitiligo have been demonstrated to exert destructive effects on melanoma cells in both in vitro and in vivo.¹⁰ Nevertheless, comprehensive research elucidating the correlation between vitiligo and melanoma at the genetic expression level remains scarce.

Transcriptomics is a technology that can resolve the dynamics of gene expression and its regulatory mechanisms. However, it is difficult to capture deeper molecular mechanisms from a single omics. Therefore, integrating RNA-seq, scRNA-seq and other multi-omics data for comprehensive analysis is the current trend.¹¹ Meanwhile, in the face of the resulting large amount of high-dimensional, heterogeneous, and non-linearly correlated data, traditional analysis methods have limited effectiveness. Machine learning algorithms can effectively mine key driving features and construct predictive models from complex multi-omics data, providing precise targets for the development of new diagnostic and therapeutic strategies.¹²

In this study, a range of bioinformatics techniques are employed to identify prognostic-related biomarkers from metastatic melanoma and vitiligo. These biomarkers are then used to construct a prognostic model, a nomogram and predict potential drugs. This study may facilitate the acquisition of novel insights into the clinical management and pathogenesis of melanoma.

Methods

Datasets of Melanoma and Vitiligo

In this study, the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) and The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/) were utilized to obtain gene expression data and clinical information on patients with malignant melanoma and vitiligo. The GSE75819¹³ and GSE46517¹⁴ datasets was employed to identify differentially expressed genes (DEGs) and to construct a gene co-expression network associated with vitiligo and metastatic melanoma, respectively. The TCGA-SKCM cohort was used to construct a prognostic model for melanoma patients. And the GSE65904¹⁵ dataset was used for the external validation of the prognostic model.

Data Processing

The datasets GSE75819, GSE46517, and GSE65904 consist of microarray data. The annotation file provided by the corresponding platform was employed to convert the probe ID to the gene symbol. If multiple probes correspond to the same gene, calculated the average of these probes. The normalizeBetweenArrays() function from the “limma”¹⁶ R package was used to perform quantile normalization while log2 scaling the data to reduce the order of magnitude difference. The TCGA-SKCM cohort was bulk RNAseq data (TPM values) scaled by log2(TPM + 1). Samples with overall survival <30 days or missing survival information were removed, resulting in the inclusion of 450 melanoma samples for the construction of the prognostic model.

Differential Expression Analysis

Differential expression analysis was performed using the “limma” R package to obtain DEGs of vitiligo relative to non-lesional skin and metastatic melanoma relative to primary melanoma from GSE75819 and GSE46517. Screening thresholds were false discovery rate (FDR) < 0.05 and |log2 fold change| > 0.585. Volcano and heatmap plots were visualized using the “ggplot2” and “pheatmap” R packages.

Weighted Gene Co-Expression Network Analysis (WGCNA)

To identify the modules and key genes associated with the disease phenotype, WGCNA was performed on GSE75819 and GSE46517 using the “WGCNA” R package. According to the scale-free topology criterion, the best softPowers were obtained using the pickSoftThreshold() function. Using the blockwiseModules() function, a co-expression network was constructed with the minimum module size set to 300 to obtain modules that were closely related to the disease phenotype. Finally, the genes were filtered from the module based on the criteria of module membership (MM) > 0.8 and gene significance (GS) > 0.5.

Vitiligo and Metastasis Related Genes (VMRGs) Identification and Functional Enrichment

The VMRGs were derived from the intersection of four gene sets: (1) DEGs obtained from GSE75819; (2) DEGs obtained from GSE46517; (3) vitiligo phenotype genes of GSE75819 obtained by WGCNA; and (4) metastatic melanoma phenotype genes of GSE46517 obtained by WGCNA. Meanwhile, the results of (1) (3) and (2) (4) were taken for intersection to obtain VMRGs. The Venn diagram of gene intersections was visualized using “ggvenn” R package. “clusterProfiler”¹⁷ R package was used to perform functional enrichment analysis (gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG)) on the VMRGs to obtain the major biological terms.

Screening Biomarkers and Building a Prognostic Model with Machine Learning

The entire TCGA-SKCM cohort was used as the training set, while GSE65904 was used as the external validation set. A univariate Cox regression was initially performed to identify VMRGs that were closely associated with overall survival via the “survival” R package. To minimize the potential bias in target screening, two machine learning algorithms were employed to further screen hub VMRGs. Based on the results of the univariate Cox regression, the least absolute shrinkage and selection operator (LASSO) Cox regression was used to further screen genes associated with prognosis via the “glmnet” R package. Additionally, the random forest (RF) algorithm was also employed to further narrow down the range of prognostically relevant genes via the “randomForest” R package. The hub VMRGs that were identified as significant in the LASSO-Cox and RF screening procedures were incorporated into the construction of the prognostic model. Multivariate Cox regression was employed to ascertain the independent prognostic factors and to develop the prognostic model via the “survival” R package. The risk scores were calculated using the prognostic model, and the optimal cutoff value was determined using the survminer R package to differentiate between high-risk and low-risk groups. The following formula was employed for the calculation of risk scores:

The disparity in survival prospects between the various risk groups was evaluated through the application of the Kaplan-Meier (KM) method. The “ggrisk” R package was used to visualize the distribution of risk scores and survival status. The receiver operating characteristic (ROC) curve was plotted for 1-, 3-, and 5-year time points using the “timeROC” R package, and the area under the curve (AUC) was calculated.

Establishment and Validation of the Nomogram

Univariate and multivariate Cox regression were used to ascertain whether risk scores were independent prognostic factors in comparison to age, sex, T stage, N stage, and M stage. Meanwhile, a nomogram was constructed based on risk scores, age, sex, T stage, N stage, and M stage to predict 1-, 3-, and 5-year time points overall survival of melanoma patients by the “rms” R package. Calibration curve analysis was employed to evaluate the accuracy and calibration of the nomogram based on the KM method. To assess the clinical value of the nomogram, decision curve analysis (DCA) was conducted via the “ggDCA” R package.

Gene Set Enrichment Analysis (GSEA)

GSEA¹⁸ is a commonly employed method for determining whether a predefined set of genes exhibits a statistically significant and consistent difference between two biological states. The differential expression analysis between the high-risk and low-risk groups was performed based on FDR < 0.05 and |log2 fold change| > 0.585 via the “limma” R package. The DEGs were subsequently employed for GSEA. The Hallmark gene sets (h.all.v2024.1.Hs.symbols) from the Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb) were utilized as reference gene sets. GSEA was conducted using the “fgsea” R package, and meaningful enriched pathways were selected based on adjusted p-value < 0.05 and sorted according to normalized enrichment score (NES).

Assessment of the Immunological Properties in Melanoma

The level of infiltration of 22 different immune cells between high-risk and low-risk groups was assessed using the CIBERSORT¹⁹ method. ESTIMATE was employed to calculate the stromal scores and immune scores between the two groups. The difference in the expression of immune checkpoint genes (ICGs) was also evaluated between the two groups. Pearson correlation analysis was utilized to calculate risk scores as well as the correlation of hub VMRGs with immune cells and ICGs.

Single-Cell RNA Sequencing (scRNA-Seq) Analysis

The raw data were obtained from three cutaneous melanoma samples from the GSE215120²⁰ dataset, and the “Seurat” R package was employed for scRNA-seq analysis. To exclude low-quality cells, overlapping multiple cells and dead, senescent cells, cells that met the following criteria were filtered out: the number of genes detected per cell exceeded 5,000 or was below 200, the unique molecular identifier (UMI) was above 20,000, and UMI from the mitochondrial genome was above 15%. The top 2,000 highly variable genes (HVGs) were selected for inclusion in a principal component analysis (PCA) for dimensionality reduction. The optimal number of principal components was selected based on the variance explained plot. Since the data came from three different samples, the “harmony” R package was used to remove batch effects. Single-cell clustering visualization was performed using 2D uniform manifold approximation and projection (UMAP). The “clustree” R package was employed for the purpose of computing the optimal resolution in the UMAP reduction process. Each cluster then was annotated with the well-known cell marker genes. Subsequently, the expression of hub VMRGs was investigated in melanoma single-cell samples. The “irGSEA” R package was employed to calculate the AUCell scores of hub VMRGs. The “CellChat” R package was employed to further elucidate the communication relationships between cells.

Small-Molecule Compounds Prediction and Molecular Docking

The Connectivity Map (CMap) database (https://clue.io/) was used to predict small-molecule compounds with therapeutic effects based on DEGs between high risk and low risk groups. The structures of compounds and proteins were obtained from PubChem (https://pubchem.ncbi.nlm.nih.gov/) and RCSB PDB (https://www.rcsb.org/) databases, respectively. Preprocessing of proteins such as remove ligands, waters, etc. and visualization of results were performed using PyMol Open-Source v3.0. Subsequently, molecular docking was performed using Dockey v1.0.3²¹ to verify the binding performance of compounds and proteins. A binding affinity < −5.0 kcal/mol indicated that the small molecule ligand exhibits favourable binding activity with the receptor protein.

Statistical Analysis

All analysis was conducted using R software v4.4.1. As previously stated, all R packages employed, and the corresponding statistical methodologies were outlined in the preceding section. A p-value of less than 0.05 was deemed to indicate a statistically significant difference.

Results

Vitiligo-Related Features and Functions

The GSE75819 dataset was normalized and processed before differential expression analysis. The sample distributions were visualized using boxplots (Supplementary Figure S1A and B) and PCA plots (Supplementary Figure S1C and D). The differential expression analysis of the processed samples resulted in 1,881 up-regulated genes and 1,800 down-regulated genes (FDR < 0.05 and |log2 fold change| > 0.585). The volcano plot and heatmap were employed to visualize these DEGs (Figure 1A and B). The WGCNA, based on the expression data of GSE75819, was performed to identify the key gene modules and genes associated with the vitiligo phenotype. Following clustering of the samples, outliers were identified, and 100 were used as the threshold for their removal. This resulted in the removal of a total of 5 samples (Supplementary Figure S1E and F). After the screening process, softPower 3 was identified as the optimal value for the network (Figure 1C). A total of 10 modules were ultimately clustered in GSE75819 (Figure 1D and E), among which the turquoise module exhibited a high correlation with the vitiligo phenotype (−0.92, p < 0.05; Figure 1E). Therefore, the turquoise module was selected for further analysis. Genes that satisfied the criteria of MM > 0.8 and GS > 0.5 were considered statistically significant (Figure 1F), resulting in the identification of 1,819 co-expressed genes associated with the vitiligo phenotype. The intersection of the DEGs and the co-expressed genes was taken, and 1,555 vitiligo-related genes were obtained (Figure 1G). These genes were analyzed by GO and KEGG. The biological process (BP) was primarily associated with mitochondrial translation, whereas the cellular component (CC) was linked to mitochondrial protein-containing complex. And the molecular function (MF) was associated with the structural constituent of ribosome (Figure 1H). The KEGG pathway of genes associated with vitiligo was found to be enriched for pathways related to the ribosome and oxidative phosphorylation (Figure 1I).

Figure 1 Vitiligo-related features and functions. (A) Volcano plot of differential expression analysis of GSE75819. (B) Heatmap of differential expression analysis of GSE75819. (C) WGCNA screening for optimal softPower in GSE75819. (D) Dendrogram of co-expression network after merging modules in GSE75819. (E) Heatmap of correlation between modules and phenotypes. (F) Scatterplot for the screening of genes that satisfy MM > 0.8 and GS > 0.5 within the turquoise module. (G) Venn plot of DEGs and co-expressed genes of vitiligo with an intersection. (H) Barplot of GO enrichment analysis of vitiligo-related genes. (I) Barplot of KEGG enrichment analysis of vitiligo-related genes.

Metastatic Melanoma-Related Features and Functions

Similarly, the GSE46517 dataset processing results were shown in the boxplots (Supplementary Figure S2A and B) and PCA plots (Supplementary Figure S2C and D). Metastatic melanoma had a total of 1,354 DEGs (FDR < 0.05 and |log2 fold change| > 0.585) compared to primary melanoma, with 569 up-regulated and 785 down-regulated. Volcano plot and heatmap were also used to visualize the DEGs (Figure 2A and B). A WGCNA was conducted on GSE46517 to identify genes associated with metastatic melanoma. Twenty-two outlier samples were excluded based on a threshold of 110 (Supplementary Figure S2E and F), after which 5 was selected as softPower to construct the network (Figure 2C). A total of 7 modules were clustered in GSE46517 (Figure 2D and E), with the brown module being the most relevant to metastatic melanoma (−0.86, p < 0.05; Figure 2E). A total of 158 genes in the brown module satisfied the criteria of MM > 0.8 and GS > 0.5 (Figure 2F) and were consequently considered to be co-expressed genes associated with the metastatic melanoma phenotype. Subsequently, a total of 156 metastatic melanoma-associated genes (Figure 2G) were obtained after intersection of DEGs and co-expressed genes, which were analyzed by GO and KEGG. GO analyses were primarily concerned with the epidermis and skin development (Figure 2H). The KEGG analysis definitively showed that metastatic melanoma-associated genes were enriched for only two pathways: Staphylococcus aureus infection and estrogen signaling pathway (Figure 2I).

Figure 2 Metastatic melanoma-related features and functions. (A) Volcano plot of differential expression analysis of GSE46517. (B) Heatmap of differential expression analysis of GSE46517. (C) WGCNA screening for optimal softPower in GSE46517. (D) Dendrogram of co-expression network after merging modules in GSE46517. (E) Heatmap of correlation between modules and phenotypes. (F) Scatterplot for the screening of genes that satisfy MM > 0.8 and GS > 0.5 within the brown module. (G) Venn plot of DEGs and co-expressed genes of metastatic melanoma with an intersection. (H) Barplot of GO enrichment analysis of melanoma metastasis-related genes. (I) Barplot of KEGG enrichment analysis of melanoma metastasis-related genes.

Machine Learning Hub VMRGs Identification and Prognostic Model Construction

Following the intersection of genes associated with both vitiligo and metastatic melanoma, a total of 27 VMRGs were identified (Figure 3A). To obtain prognostically relevant VMRGs, we conducted a further screening of 18 VMRGs that were significantly associated with overall survival in melanoma by univariate Cox regression in the TCGA-SKCM cohort (HR > 1, p < 0.05; Figure 3B). Subsequently, the genes were screened with two classical machine learning algorithms to enhance reliability and accuracy. Based on the LASSO-Cox regression results, 12 genes were selected (Figure 3C and D). In addition, MeanDecreaseAccuracy (MDA) scores in RF revealed the importance of each gene, and we selected the top10 MDA genes (Figure 3E). The LASSO Cox regression and RF overlapped resulted in a total of 6 hub VMRGs (DUOX1, GJB3, NOTCH3, PKP1, PTK6, and PTPRF). Then, we constructed a prognostic model based on 6 hub VMRGs using multivariate Cox regression. The training set was derived from 450 samples from the TCGA-SKCM cohort, while the external validation set was derived from 214 samples from GSE65904. The risk score formula for the model was as follows:

Figure 3 Construction of prognostic models by machine learning. (A) Venn plot of DEGs and co-expressed genes from vitiligo and metastatic melanoma, respectively with an intersection. (B) Univariate Cox regression forest plot. (C and D) The LASSO-COX method. (E) The RF algorithm with genes ranked based on the MDA and MDG. (F) Multivariate Cox regression forest plot.

Among the 6 hub VMRGs, GJB3, NOTCH3, PKP1, and PTK6 were independently associated with prognosis (Figure 3F), with GJB3 being a protective factor (HR < 1, p < 0.05) and NOTCH3, PKP1, and PTK6 being risk factors (HR > 1, p < 0.05).

Validation of the Prognostic Model

The training set was divided into two groups, designated as high-risk (n = 299) and low-risk (n = 151), based on the calculated optimal cutoff value of 0.796. Kaplan-Meier analysis revealed that patients in the high-risk group exhibited a significantly poorer overall survival (p < 0.05; Figure 4A). The AUC values of the 1-, 3-, and 5-year time points ROC curves were 0.718, 0.667, and 0.646, respectively. These values indicate that the prognostic model maintained good sensitivity and specificity over long latitudes (Figure 4B). To assess the reliability of the prognostic model, risk scores for the external validation set were calculated using the same risk score formula. In the external validation set, the high-risk (n = 160) and low-risk (n = 50) groups were divided based on the optimal cutoff value of 2.93. Kaplan-Meier analysis also demonstrated that the high-risk group exhibited a poorer prognosis (p < 0.05; Figure 4D). The 1-, 3-, and 5-year time points ROC (AUC) for the external validation set were 0.617, 0.605, and 0.596, respectively (Figure 4E). The distributions of risk score, survival status, and expression levels of the hub VMRGs for the training set and external validation set were illustrated in Figure 4C and F, respectively. Ultimately, the subgroups were delineated according to age, sex, T, N, and M stage in the TCGA-SKCM cohort. Kaplan-Meier analyses revealed that the high-risk group had a lower overall survival rate than the low-risk group, except for the M1 subgroup (Supplementary Figure S3A–J). It’s worth noting that GSE65904 includes more advanced-stage patients, potentially limiting model generalizability. It is of note that the model demonstrated prognostic capability, indicating potential use in predicting melanoma progression in patients.

Figure 4 Validation of the prognostic model. (A) Kaplan-Meier plots in the training set (TCGA-SCKM). (B) ROC curves in the training set (TCGA-SCKM). (C) The distribution of risk score, associated survival data, and expression of the hub VMRGs in the training set (TCGA-SCKM). (D) Kaplan-Meier plots in the external validation set (GSE65904). (E) ROC curves in the external validation set (GSE65904). (F) The distribution of risk score, associated survival data, and expression of the hub VMRGs in the external validation set (GSE65904).

Nomogram Construction and Validation Based on Risk Score

To ascertain whether the risk score is an independent prognostic factor for melanoma patients, univariate and multivariate Cox regression were performed on the TCGA-SKCM cohort, respectively. The univariate Cox regression (Figure 5A) revealed that the risk score, age, T4, N1, N3 and M were correlated with prognosis (HR > 1, p < 0.05). Meanwhile, in the multivariate Cox regression (Figure 5B), the risk score, age, T4, N1, and N3 were identified as independent risk factors (HR > 1, p < 0.05). By integrating the risk score with all clinicopathological features, we developed a nomogram to facilitate the optimal utilization of the risk score in a clinical setting. Figure 5C presented a nomogram and a specific example of use. The calibration plot demonstrated that the predictive outcomes of the nomogram were in close alignment with the ideal model (Figure 5D). Furthermore, the ROC curve substantiated the favorable predictive capacity of the nomogram at the 1-, 3-, and 5-year time points (Figure 5E). The DCA illustrated that the nomogram exhibited superior clinical utility at the 3- and 5-year time points (Figure 5F). The aforementioned results indicated that the nomogram developed on the basis of risk score exhibited superior predictive performance in relation to the prognosis of melanoma patients.

Figure 5 Nomogram construction and validation. (A) Univariate Cox regression forest plot. (B) Multivariate Cox regression forest plot. (C) A nomogram model was employed to predict 1-, 3-, and 5-year time points survival probabilities for melanoma patients. The red text represented an example sample’s points on each item as well as survival probabilities at different time points. (D) Calibration curves. (E) ROC curves. (F) Decision curve analysis plot. ***p < 0.0001.

Functional Annotations for Hub VMRGs Based on Risk Score

To understand the survival differences between the high-risk and low-risk groups, we identified a total of 583 DEGs (Figure 6A). Of these, 491 were upregulated and 92 were downregulated (FDR < 0.05 and |log2 fold change| > 0.585). Subsequently, the identified DEGs were subjected to GO and KEGG enrichment analyses. The GO analyses indicated that these DEGs were predominantly associated with skin and epidermal development and differentiation (Figure 6B). KEGG was primarily associated with pathways such as the cytoskeleton in muscle cells and ECM-receptor interaction, as illustrated in Figure 6C. Subsequently, GSEA based on Hallmark gene sets (Figure 6D) demonstrated that epithelial mesenchymal transition (EMT), angiogenesis, and Notch signaling pathways were significantly up-regulated, while IL6/JAK/STAT3 signaling and interferon alpha response were down-regulated. These findings may contribute to the poor prognosis observed in the high-risk group relative to the low-risk group.

Figure 6 Functional annotations for hub VMRGs based on risk score. (A) Volcano plot of differential expression analysis in the TCGA-SKCM cohort. (B) Barplot of GO enrichment analysis of genes associated with the risk score. (C) Barplot of KEGG enrichment analysis of genes associated with risk scores. (D) GSEA of DEGs in high- and low-risk groups.

The Landscape of the Immunological Properties in Melanoma

The infiltration of immune cells exerts a significant influence on tumor development. The outcomes of the CIBERSORT algorithm, as illustrated in Figure 7A and B, indicated that patients in the high-risk group exhibited elevated proportions of T cells CD4 memory resting, macrophages M0 and M2, while the proportions of plasma cells, T cells CD8, T cells CD4 memory activated, T cells gamma delta, and macrophages M1 were diminished, in comparison to the low-risk group. The ESTIMATE results indicated that the high-risk group exhibited lower immune scores and no difference in stromal scores (Figure 7C and D). This indicates that patients in the high-risk group may exhibit a suppressive cellular immune response. The vast majority of ICGs were highly expressed in low-risk groups, indicating that immune checkpoint blockade (ICB) is more beneficial (Figure 7E). Pearson correlation analyses were conducted for the hub VMRGs, risk score, immune cells, and ICGs (Figure 7F and G). The findings revealed that the NOTCH3 expression pattern in immune cells was more aligned with the high-risk group. The other hub VMRGs and risk score showed a negative correlation with most ICGs.

Figure 7 The landscape of the immunological properties in melanoma. (A) Rainbow plot of the relative percentages of each immune cell. (B) Boxplot illustrates the relative distribution of immune cells between the high- and low-risk groups. (C) Boxplot illustrates the comparison of ImmuneScore between the high- and low-risk groups. (D) Boxplot of StromalScore comparisons between the high- and low-risk groups. (E) Boxplot illustrates the comparison of ICG expression levels between the high- and low-risk groups. (F) Correlation between the hub VMRGs, the risk score and the immune cells. (G) Correlation between the hub VMRGs, the risk score and the ICGs expression levels. *p < 0.05, **p < 0.001, ***p < 0.0001, ns p > 0.05.

Expression Characteristics of Hub VMRGs and Association with Prognostic Features

The expression characteristics of the hub VMRGs in the GSE75819, GSE46517 and TCGA-SKCM cohort were further investigated (Figure 8A–C). The results indicate that hub VMRGs were less expressed in patients with vitiligo and post-metastatic melanoma, exhibiting a consistent pattern. Conversely, they were more highly expressed in the high-risk group, suggesting a worse prognosis. Kaplan-Meier plots (Figure 8D–I) indicated that high expression of all hub VMRGs except for GJB3 predicted a worse prognosis. Furthermore, there was a significant association between overall survival and NOTCH3, PKP1, and PTK6 (p < 0.05). Finally, we additionally validated in the GSE8401²² (52 metastatic melanoma samples vs 31 primary melanoma samples, Figure 8J) and GSE7553²³ (40 metastatic melanoma samples vs 14 primary melanoma samples, Figure 8K) datasets for expression levels of hub VMRGs, and the results and trends were consistent with our previous findings, suggesting that our results were reliable and stable.

Figure 8 Expression characteristics of hub VMRGs and association with prognostic features. (A) The expression levels of hub VMRGs in GSE75819. (B) The expression levels of hub VMRGs in GSE46517. (C) The expression levels of the hub VMRGs in the TCGA-SKCM cohort. (D–I) Kaplan-Meier plots demonstrate the correlation between the expression levels of GJB3, NOTCH3, PKP1, PTK6, and PTPRF in TCGA-SCKM cohort and the probability of patient survival. (J) The expression levels of hub VMRGs in GSE8401. (K) The expression levels of hub VMRGs in GSE7553. *p < 0.05, **p < 0.001, ***p < 0.0001.

Validation of Cell Specificity of Hub VMRGs by scRNA-Seq

Following quality control of the GSE215120 dataset, a total of 23,882 single cells originating from cutaneous melanoma samples were obtained (Supplementary Figure S4A and B). The top 2000 HVGs were selected for inclusion in the PCA (Supplementary Figure S4C), and the top 20 principal components were selected for UMAP reduction (Supplementary Figure S4D). Supplementary Figure S4E showed UMAP plots before and after batch effect removal. Following an evaluation of alternative resolution values, 0.5 was selected as the optimal value for clustering, resulting in the identification of 16 clusters (Supplementary Figure S4F and Figure 9A). The clusters were ultimately annotated in accordance with well-known cell marker genes (Figure 9A and B). The AUCell scores of hub VMRGs were markedly concentrated on fibroblast cells (Figure 9C and D). The analysis of cell communication revealed that the fibroblast cells exhibited the strongest outgoing interaction strength (Figure 9E and F). To this end, we conducted further investigation into the interaction between fibroblast cells and melanoma cells. Our findings revealed a significant interaction between the two cells, with COL1A1/COL1A2/COL6A1/COL6A2-CD44 being the most pronounced (Figure 9G). Finally, we investigated the correlation between the hub VMRGs and these genes (Figure 9H). The Pearson correlation demonstrated that NOTCH3 exhibited a significant positive correlation with COL1A1, COL1A2, COL6A1, and COL6A2 (p < 0.05).

Figure 9 scRNA-seq analysis. (A) UMAP reduction plot, from left to right: colored by clusters; cell types. (B) Heatmap of cell marker genes expression in different clusters and association with cell type. (C) AUCell score of hub VMRGs in the UMAP reduction plot. (D) The distribution of the AUCell score of hub VMRGs. (E) Circle plot of signaling networks in melanoma. (F) Dot plot was presented to demonstrate the strength of cellular output and input communication. (G) Cell-cell interaction between fibroblast and melanoma. The red text and boxes highlight the strong communication strength of these ligand-receptor pairs. (H) Correlation between hub VMRGs and key cell-cell interaction genes. *p < 0.05, **p < 0.001, ***p < 0.0001.

Prediction of Potential Therapeutic Compounds

150 up-regulated and 92 down-regulated genes between high risk and low risk groups were used to predict small-molecule compounds in CMap database. We selected the five compounds with the lowest scores as the most promising therapeutic molecules: methyl-angolensate, byssochlamic-acid, homoharringtonine, piperacillin and cephaeline (Figure 10A). The molecular docking results were shown in Figure 10B, where the lower the affinity (kcal/mol), the tighter the binding was considered to be. Unfortunately, the structure of GJB3 was not obtained. We then show specific results of NOTCH3 docking to byssochlamic-acid, methyl-angolensate, and cephaeline (Figure 10C–E).

Figure 10 Small-molecule compounds prediction and molecular docking. (A) Predicting results for potential therapeutic compounds by CMap database. (B) Heatmap of the affinity of molecular docking. (C–E) Detailed docking diagram of NOTCH3 docking with byssochlamic-acid, methyl-angolensate, and cephaeline, respectively.

Discussion

Melanoma is a highly malignant tumor with a high propensity for invasion and metastasis, which collectively portend a very poor prognosis.²⁴ Despite the significant advancements that have been made in the treatment of melanoma, numerous challenges remain.²⁵ There have been numerous studies showing that patients with vitiligo have a lower risk of developing melanoma and that patients with melanoma who develop vitiligo-like skin changes have a better prognosis.^5–7 Although the fact that the potential relationship between vitiligo and melanoma has now been elucidated in terms of humoral and cellular immunity,⁹ there remains a paucity of research investigating the crosstalk between vitiligo and melanoma at the gene expression level, as well as studies based on vitiligo-related biomarkers.

In this study, we screened genes associated with vitiligo and metastatic melanoma combined with differential expression analysis and WGCNA. Furthermore, hub VMRGs (DUOX1, GJB3, NOTCH3, PKP1, PTK6, and PTPRF) were identified as being significantly associated with prognosis by machine learning algorithms, and validated the expression levels in independent cohorts. A prognostic model based on these genes in the TCGA-SKCM cohort was developed and also externally validated in GSE65904. The prognostic model classifies melanoma patients into high-risk and low-risk groups. The high-risk group had a worse prognosis. Consequently, A nomogram was developed to provide new insights for clinical decision-making. Finally, the immune landscape, single-cell characteristics and therapeutic compounds of melanoma patients were subjected to further investigation. Compared with previous studies,^26,27 this study is the first to use multi-omics and machine learning to integrate vitiligo data and discover Notch signaling pathway unique crosstalk mechanism.

Among the six hub VMRGs, firstly, NOTCH3 and Notch signaling pathway may play a pivotal role in the crosstalk between vitiligo and metastatic melanoma. Vitiligo is a condition characterized by the loss of pigmentation in the skin, caused by the dysfunction or absence of melanocytes. The development and maintenance of the melanocyte lineage is dependent on Notch signaling.²⁸ Inactivation of Notch signaling eliminates melanocytes and melanocyte stem cells, resulting in severe hair pigmentation defects.²⁹ The inactivation of NOTCH1 in vitiligo has been shown to result in the loss of epidermal and/or follicular active melanocytes³⁰ Despite the absence of prior reports associating NOTCH3 with vitiligo, our study revealed a notable reduction in NOTCH3 expression in vitiligo tissues relative to normal skin tissues (Figure 8A). Conversely, aberrant activation of the Notch signaling pathway has been observed to promote melanoma development and metastasis. Furthermore, inhibition of the Notch signaling pathway has been demonstrated to suppress melanoma cell proliferation.³¹ NOTCH3 has been shown to mediate communication between melanoma cells and endothelial cells, as well as tumor cell migration.³² In our study, NOTCH3 was found to be highly expressed in the high-risk group (Figure 8C), indicating a worse prognosis (Figure 8F). Furthermore, there was a strong correlation between NOTCH3 expression and immunosuppression (Figure 7F and G). GSEA demonstrated that the Notch signaling pathway was up regulated in the high-risk group (Figure 6D). In conclusion, there is evidence to suggest that NOTCH3 plays a pivotal role in the pathogenesis of both vitiligo and melanoma.

Additionally, other hub VMRGs are also linked to the development and progression of cancers. PKP1 has significant effects on melanoma cell proliferation, migration, invasion and cell cycle³³ PTK6 has been identified as a pro-carcinogenic factor in a number of cancer types, including breast cancer³⁴ and prostate cancer.³⁵ Additionally, there is evidence to suggest that it may play a role in the development of metastatic melanoma.³⁶ GJB3 has been demonstrated to promote liver metastasis in pancreatic cancer.³⁷ The potential of this marker to be a new melanoma marker was first reported in 2019, but further studies are currently lacking.³⁸ It was observed that DUOX1 exerts a dual role in melanoma, inhibiting proliferation but promoting metastasis.³⁹ PTPRF has manifested inhibitory effects against gastric adenocarcinoma⁴⁰ and breast cancer,⁴¹ yet the mechanism of action in melanoma remains opaque.

The infiltration of immune cells represents a pivotal mechanism in the progression of tumors. The proportion of plasma cells, CD8 T cells, memory-activated CD4 T cells, and gamma delta T cells infiltrated in the high-risk group of patients was found to be significantly reduced (Figure 7B), indicating that T-cell immunity were suppressed in the high-risk group. We observed a rise in M2 macrophage infiltration and a lower degree of M1 macrophage infiltration in the high-risk group (Figure 7B). These results suggest that tumors in the high-risk group have a greater tumor-associated macrophage (TAM) infiltration. ICB has played a significant role in the treatment of melanoma in recent years. A comparison of the expression levels of ICGs between high- and low-risk groups revealed that the majority of ICGs were relatively highly expressed in the low-risk group (Figure 7E), indicating that the low-risk group may potentially derive greater benefit from ICB.

In recent years, the development of scRNA-seq techniques has facilitated a deeper understanding of cellular composition and interactions within tumors. CD44 is a critical mediator in the EMT process observed in tumor development.⁴² Our findings indicate that hub VMRGs may play a role in melanoma proliferation and metastasis by influencing fibroblasts to promote EMT (Figure 9C–G). The results of the correlation analysis indicated a significant and positive correlation between NOTCH3 and COL1A1, COL1A2, COL6A1, and COL6A2, which are collagen genes associated with fibroblasts.⁴³

Finally, we found that methyl-angolensate, byssochlamic-acid, homoharringtonine, piperacillin and cephaeline may be potential therapeutic compounds for targeting hub VMRGs based on Cmap database and molecular docking. In these compounds, homoharringtonine can inhibit melanoma cell proliferation by inducing DNA damage, apoptosis and G2/M cell cycle arrest.⁴⁴ Meanwhile, it has been indicated to inhibit the NOTCH/MYC pathway,⁴⁵ which is consistent with our results. Although the remaining compounds are not reported to be directly associated with melanoma, they are closely related to tumor therapy and require further research in the future. Taken together, our results suggest that NOTCH3 may be a new potential therapeutic target for the treatment of melanoma.

It should be noted that our study is not without limitations. We have developed a prognostic model using publicly available data. While external validation indicates a degree of robustness, further clinical validation is required to confirm its generalizability and to address the possibility of overfitting to specific sample sets.

Conclusions

We firstly investigated the crosstalk mechanisms between vitiligo and metastatic melanoma by integrating multi-omics and machine learning. Prognostic model and nomogram were constructed from hub genes (DUOX1, GJB3, NOTCH3, PKP1, PTK6 and PTPRF), which was able to accurately stratify risk for melanoma patients. NOTCH3 and Notch signaling pathway may play a pivotal role in the crosstalk between vitiligo and metastatic melanoma, but extensive human and animal experiments are needed for further validation in the future. Finally, we predicted potential small-molecule compounds and verified the affinity of hub VMRGs by molecular docking. In conclusion, this study may provide important information to improve the prognosis of melanoma patients.