3.1 Genomic mapping differences between normal and CM
The study’s flowchart is depicted in Supplementary Fig. 1. Using cluster analysis and principal component analysis (PCA), we first demonstrated the genetic differences between normal skin tissues and CM (Fig. 1A-B). By comparing 469 CM and 556 normal skin tissues, we discovered that the expression of 4555 genes was significantly altered in tumor tissues against normal tissues (P < 0.01, |logFC|≥2), with 2296 genes considerably up-regulated and 2259 genes significantly down-regulated (Fig. 1C). The GO enrichment analysis demonstrates that CM genetic variants play a role in TME immune components and matrix-related biological processes, such as cytokine-mediated signaling pathway, cell chemotaxis, and positive regulation of response to external stimulation (Fig. 1D-G). The KEGG pathway indicates that these genes are also involved in immune-related signaling pathways, such as cytokine-cytokine receptor interaction, viral protein interaction with the cytokine-cytokine receptor, viral protein interaction with cytokine and cytokine receptor, and natural killer cell-mediated cytotoxicity (Fig. 1H). The heatmap results also imply that immunological and matrix-related pathways play a significant role in the CM genome (Fig. 1I). The PPI network reveals a close interaction between CM-related genes at the protein level (Fig. 2A). We identified a sub-network module with solid predictive value, in which CD86, CXCL9, FCGR3A, GZMB, PRF1, STAT1, and TLR7 were designated as key molecules with high connectedness. We also investigated the protein correlation between these key molecules (Fig. 2B-C). We found that the expression of critical molecules is significantly elevated in CM samples relative to normal samples (Fig. 2D). We performed dimensionality reduction using PCA to determine if these critical molecules can distinguish CM samples from normal samples. We discovered two completely disjoint populations, indicating that the expression patterns of critical molecules in normal and CM samples are distinct (Fig. 2E). In addition, we demonstrated the interaction of numerous immune molecules in the CM TME, the signaling cascade and transmission, and the distinct regulation patterns between molecules (Fig. 2F).
3.2 Multi-omics analysis to identify key molecules and mutation and survival analysis
We obtained immunohistochemistry results for seven critical molecules from the HPA database and qualitatively found protein-level expression variations between normal tissues and CM samples (Fig. 3A-B). In a subsequent confirmation, we evaluated the expression of these important molecules in cell lines and 20 pairs of tumors and surrounding normal tissues to validate the differential and significant expression of essential molecules in CM tissues. Spearman correlation analysis identified a significant positive correlation between important molecules and a strong interaction between these molecules (Fig. 3C). The mutations of CM’s essential molecules were then studied. Among the 467 patients with complete clinical annotation information, 72 patients (15.42%) had mutations in at least one gene locus. CD86 had the highest frequency of mutations in CM samples, followed by TLR7, and all essential molecules were found with gene mutation (Fig. 3D). The expression and mutation of these essential molecules may play a crucial role in the growth and metastasis of CM, as deduced by our findings. In addition, we probed into the predictive value of key molecules based on an independent CM cohort from the TCGA database using survival analysis. The TCGA-CM cohort also demonstrated variations in the expression of essential molecules between normal and tumor samples. For survival analysis, 467 CM patients with complete clinical annotation were available. Patients with high expression of CD86, CXCL9, FCGR3A, GZMB, PRF1, STAT1, and TLR7 had a significant survival benefit over those with low expression (Fig. 4A-G). Seven critical molecules displayed significantly greater expression levels in CM, and all seven key molecules as protective molecules significantly increased CM patient survival.
3.3 Evaluation of immune cell infiltration characteristics of tumor microenvironment
To further investigate the role of identified critical molecules in TME immune cell infiltration in CM patients, we analyzed the infiltration of 28 types of TME cells in normal and tumor tissues (Fig. 5A). T helper cells (type 1 and 2), activated B cells, CD4+ T cells, CD8+ T cells, immature B cells, regulatory T cells, natural killer cell, activated dendritic cell, plasmacytoid dendritic cell, MDSC, monocyte, memory B cell, macrophage, gamma delta T cell, effector memory CD4+ T cells, CD56 dim natural killer cell, immature dendritic cell, eosinophil, CD8+ T cells were highly plentiful in tumor tissue. However, other cell subsets were notably abundant in normal tissue. Then, using the PCA algorithm, we compared the infiltration patterns of TME cells in normal and tumor tissues to determine if there were any differences. After dimensionality reduction, the results demonstrated the existence of two distinct populations of TME cells (Fig. 5B). Using the ESTIMATE algorithm, we determined the immunological and mesenchymal activity in the CM microenvironment. It was discovered that immune and mesenchymal activities were much higher in tumor tissues than normal skin tissues (Fig. 5C-D). To examine the link between critical molecules and immune cells in the TME, we correlated key molecules with cellular fractions in the TME. Spearman correlation analysis revealed that these molecules were strongly positively linked with most TME cellular fractions, except T helper cells, CD56 bright natural killer cells, and neutrophil cell infiltration (Fig. 5E). In addition, the expression of seven essential molecules demonstrated a significant positive association with PD-L2 and PD-L1 and a negative correlation with CTLA4 (Fig. 5F). CTLA4 had the most significant connection with STAT1 and TLR7 (Fig. 5G-H).
3.4 Correlation model construction for prognosis and immunotherapy based on key molecules
We incorporated patient prognostic information and TME immune cell infiltration status to build the riskScore model, and we integrated the role of these essential molecules using LASSO Cox regression. RiskScore was determined by the expression of the four most representative important molecules, according to the findings (Fig. 6A-B). Based on the critical value of -7.07 (Fig. 6C) computed by the MaxStat R package, we classified the patients into high-risk and low-risk groups. We observed that the low-risk group had a considerable survival advantage over the high-risk group (Fig. 6D). In addition, the expression of these essential molecules is much higher in low-risk tumors than in high-risk, suggesting that these key molecules have a protective role in the low-risk group, consistent with the findings of our earlier investigation (Fig. 6E-F). With rising risk, patient mortality might climb significantly (Fig. 6G-H). Our examination of multivariate COX regression models incorporating basic clinical and pathological information about the patients demonstrated that riskScore could be an independent and robust predictive biomarker to evaluate CM patients (Fig. 7A). In addition, we developed a nomogram that combines the riskScore with independent clinical prognostic indicators to estimate the likelihood of patient mortality (Fig. 7B). The calibration plots demonstrated that the generated nomogram had a superior prediction ability (Fig. 7C). We displayed ROC curves based on TCGA data with AUC of 0.757, 0.675, and 0.657 for 1, 2, and 3 years, indicating the riskScore’s predictive performance is acceptable and can provide a reference for clinical decision-making (Fig. 7D). In addition, our findings demonstrated that riskScore surpassed other clinical factors, such as age and number of nodules, in predicting OS in patients with CM (Fig. 7E).
We utilized gene set enrichment analysis (GSEA) to investigate the activated biological pathways in the low-risk and high-risk groups. Compared to the low-risk group, cancer-related pathways such as P53, PI3K-AKT-mTOR, NOTCH, and WNT were considerably activated in the high-risk group (Fig. 7F-K). Then we analyzed the difference in TME cell infiltration between the low-risk and high-risk groups, and we discovered that all immune infiltrating cells, except for CD56 dim natural killer cells and CD56 bright natural killer cells, were significantly higher in the low-risk group than in the high-risk group (Fig. 8A). By correlation analysis, we discovered that riskScore values were strongly and positively linked with the majority of TME cell infiltrating rates (Fig. 8B). We also discovered a significant and positive correlation between riskScore values and the expression of immune checkpoint molecules, indicating the potential predictive role of riskScore in predicting clinical response to immunotherapy and providing a foundation for developing novel immunotherapies (Fig. 8C). As immune checkpoint blockade (ICB) has made advances in the treatment of CM over the past few years, we verified riskScore’s ability to predict the clinical response of patients to ICB therapy. In the IMvigor210 cohort treated with anti-PD-L1 therapy, low-risk patients had a considerable clinical benefit and prolonged survival (Fig. 8D). The patients with complete remission (CR) or stable disease (SD) had a lower risk (Fig. 8E). In addition, we noticed that low-risk individuals responded considerably better to PD-L1 blocking therapy than high-risk patients (Fig. 8F).
3.5 Chemotherapy drug sensitivity analysis, small molecule drug screening and molecular docking validation
We analyzed 20 common chemotherapeutic and targeted medicines and found substantial variations between the high-risk and low-risk categories in IC50 values (Supplementary Fig. 2). The results indicate that our riskScore signature can uncover prospective biomarkers of chemotherapy and targeted medication sensitivity. Then, we calculated the connection between medication-treated expression profiles and highly up-regulated expression profiles of seven key genes using the Cmap database. We then identified the top ten pharmaceuticals with negative correlations as potential treatment candidates (Table 1). Figure 9A-J shows the chemical structures of these ten compounds. AGI-6780 and Zofenopril-calcium bind well to GZMB, indicating that these two small compounds can be employed as possible target medicines to target GZMB. In addition, we utilized Pymol to generate a heatmap of the binding of CD86, FCGR3A, STAT1, TLR7, and GZMB proteins to the most strongly bound small molecules or the top two most strongly bound small molecules (Fig. 9K). The results demonstrated that the small molecules of CD86 bound to Baricitinib formed hydrogen bonds with THR-69, SER-67, GLN-16. The binding of FCGR3A to Zofenopril-calcium formed hydrogen bonds with HIS-111 and ARG-109. The binding of GZMB to Zofenopril-calcium formed hydrogen bonds with LYS-113 and ARG-87, and no hydrogen bonds were formed in the binding of AGI-6780. Small molecules in the binding of STAT1 to Calcipotriol formed hydrogen bonds with GLU-353 and GLN-271, and small molecules in the binding of TLR7 to Zofenopril-calcium formed hydrogen bonds mainly with LYS-464. The majority of receptor and ligand binding energies are less than − 7 kcal.mol− 1, indicating that the target protein and active ingredient can bind spontaneously with high affinity and stable conformation, and thus small molecule medicines are likely to act on these targets. To illustrate the molecular interactions, we chose the small molecule medication with the lowest binding energy to dock the target for docking visualization (Fig. 10A-H).
Table 1
Cmap name
|
N
|
Celline
|
Enrichment
|
FDR_Q_nlog10
|
Gabapentin
|
2
|
YAPC
|
-0.94
|
15.65
|
Baricitinib
|
3
|
HBL1
|
-0.92
|
15.65
|
DPN
|
3
|
A549
|
-0.91
|
15.65
|
AGI-6780
|
2
|
PC3
|
-0.9
|
15.65
|
Fusaric-acid
|
3
|
SKB
|
-0.9
|
15.65
|
Ru-24969
|
3
|
MCF7
|
-0.89
|
15.65
|
Calcipotriol
|
2
|
HCC515
|
-0.89
|
15.65
|
Fenoterol
|
2
|
YAPC
|
0.89
|
15.65
|
Zofenopril-calcium
|
2
|
JURKAT
|
-0.89
|
15.65
|
RS-102895
|
3
|
A549
|
-0.89
|
15.65
|
3.6 Expression validation of key molecules
We validated the differential expression of these proteins in normal skin tissues and CM using semi-quantitative analysis after selecting one appropriate field of view (the first column of Fig. 11). qRT-PCR was subsequently utilized to confirm the differential expression of these seven essential genes in cell lines and human specimens. The mRNA expression of seven key genes was considerably higher in the A375 and SK-Mel-14 cell lines than in the Hacat and PIG1 cell lines. (the second column of Fig. 11). Similarly, we discovered that the patients’ mRNA expression of seven key genes was much higher in CM than in normal skin tissue (the third column of Fig. 11). Using experimental validation at the mRNA and protein levels, we determined that CD86, CXCL9, FCGR3A, GZMB, PRF1, STAT1, and TLR7 were differentially expressed in normal skin tissues and CM and inferred that these key molecules could be potentially critical targets for the treatment of CM.