DEGs selection
After standardization and identification of the microarray results, DEGs (4640 in GSE15065,1096 in GSE46517 and 1895 in GSE114445) were selected. The overlap among three datasets included 308 genes as shown in the Venn diagram (Fig. 1a, Supplementary table 1) between primary melanoma and normal skin.
GO and KEGG Enrichment Analysis of the DEGs
Functional and pathway enrichment of DEGs was performed by DAVID, then displayed in bubble charts by R software (p<0.05, Fig. 2). GO analysis indicated that changes in biologic processes were significantly enriched in the inflammatory response, immune response, regulation of transcription. As for cellular component, DEGs were particularly enriched in extracellular region, extracellular space and extracellular exosome. Changes in molecular function were mostly enriched in chemokine activity, transcriptional activator activity, protein binding and CXCR3 chemokine receptor binding. KEGG pathway analysis demonstrated that DEGs played pivotal roles in pathways in cytokine-cytokine receptor interaction, chemokine signaling pathway and pathways in cancer. These results revealed that immune response and inflammatory response play important roles in SKCM tumorigenesis.
PPI network construction and hub genes identification
The PPI network was constructed via Cytoscape, including 247 nodes and 792 edges (Fig. 1b). Then, the most important PPI network module was obtained using MCODE, consisted of 12 nodes and 64 edges (Fig. 1c). After selection, CCL5, PTGER3, IL6, CXCL13, CCL27, CCL4, CXCL9, NMU, CXCL2, GAL, NPY1R, CXCL10 were considered as hub genes of the network.
Function analysis of hub genes
As expected, functional annotation obtained from Metascape suggested that hub genes were mainly enriched in regulation of T cell chemotaxis, peptide ligand-binding receptors and chemokine receptors bind chemokine signaling (Fig. 3a-b, P<0.001). Similarly, ClueGO (Fig. 3c, P<0.01) revealed that the most involved pathways were interleukin-10 signaling, peptide ligand-binding receptors and chemokine receptors bind chemokine signaling.
The prognostic and diagnostic value of hub genes in SKCM patients
Using the Kaplan‑Meier method, the prognostic values of the hub genes in SKCM patients were determined. Seven hub genes were significantly associated with OS as shown in Fig. 4. High expression of CXCL9 (HR= 0.6; P=2E-4), CXCL10 (HR=0.57; P=2.8E-5), CXCL13 (HR=0.56; P=2.1E-5), CCL4 (HR=0.48; P=8.2E-8), CCL5 (HR=0.53; P=3.2E-6) were significantly associated with better OS. In contrast, NMU (HR=1.6; P=3E-4) and GAL (HR=1.6; P=0.001) were found associated with poor overall survival. Furthermore, CCL4 (HR=0.76; P=0.026), CCL5 (HR=0.76; P=0.027) were also found significantly associated with better disease-free survival (Fig. 4). Except GAL (AUC=0.536), the rest of hub genes, including CCL4 (AUC=0.872), CCL5 (AUC=0.775), CXCL9 (AUC=0.834), CXCL10 (AUC=0.786), CXCL13 (AUC=0.807) and NMU (AUC=0.829), were of diagnostic value (p<0.001) (Fig. 5).
Validation the aberrant expression of prognostic genes and their clinical characteristics
Based on data from GEPIA database, all of the prognostic genes were found differential expressed in SKCM vs normal skin tissue (P<0.05, Fig. 6).
In addition, CXCL9 (P=7.74E-5), CXCL10 (P=1.05E-4), CXCL13 (P=8.15E-4), CCL4 (P=0.002), CCL5 (P=0.001) were significantly correlated to SKCM pathological stages (Fig. 7). Among them, CXCL9, CXCL10, CXCL13, CCL4, CCL5 were higher expressed in stage I and decreased in the subsequent stages.
Subsequently, all of the samples from TCGA were ordered according to the expression level of CXCL9, and statistical analysis was conducted to compare the head 80 with the tail 80 samples. The heatmap (Fig. 8b) and statistical results suggested that higher expression levels (tail 80 samples) of CXCL9, CXCL10, CXCL13, CCL4, CCL5 were related to more patients of lower Breslow depth (0-1.5mm, 51.2%), BRAF mutation (60%) and lower Clark levels (Ⅰ-Ⅲ 45% and Ⅳ-Ⅴ, 55%), reduced T4 stage (T4, 25%). With the decrease levels of these five chemokines (head 80 samples), there were more patients with higher Breslow depth (0-1.5mm, 15% and >1.5mm, 85%), BRAF WT (wild type, 56.3%), higher Clark levels (Ⅰ-Ⅲ 17.5% and Ⅳ-Ⅴ, 82.5%) and T4 stage (T4, 47.5%).
All of the data above suggested that five chemokine family members play critical roles in cutaneous melanoma tumorigenesis and progression.
Genetic alteration, correlation coefficient and key transcription factors analyses of prognostic genes in patients with SKCM
As a result, there were nearly 6% (CXCL9), 5% (CXCL10), 3% (CXCL13), 4% (CCL4), 5% (CCL5), 3% (NMU), 4% (GAL) of SKCM samples had genetic alteration. The most common genetic change among five chemokines (CXCL9, CXCL10, CXCL13, CCL4, CCL5) was enhanced mRNA expression (Fig. 8a). While genetic change in GAL and NMU were mainly related to amplification.
Then, we assessed the correlation coefficients in prognostic genes. As expect, significant positive correlations were observed between five chemokines (Fig. 8c). Among them, CXCL9 and CXCL10 had the highest positive correlation with a Spearman's correlation coefficient of 0.92. CXCL9 exhibited the tightest association with all other chemokines, with a median Spearman's correlation coefficient of 0.85.
Using TRRUST database, we explored potential transcription factor targets of the five chemokines. As a result, IRF7, IRF3 (FDR=3.09E-05); RELA, NFKB1 and IRF1 (FDR=1.5E-04) were found to be the key transcription factors for CCL4, CCL5, CXCL10 (Table 2).
Validation of five chemokines expression in SKCM tissues from the FAHSU cohort
Therefore, IHC was performed to validate the expression of five chemokine family members. Consistent with the mRNA expression, we found significantly elevated CXCL9, CXCL10, CXCL13, CCL4 and CCL5 proteins expression in the SKCM than in the normal tissues. The results and the plots of IHC score were illustrated in Fig. 9.
Immune infiltration analysis
We performed correlation analysis between prognostic genes expression and immune infiltration level for SKCM. There was a positive correlation between CXCL9 expression and the infiltration of B cells (Cor=0.238, p=3.58e−07), CD8+T cells (Cor=0.643, p=1.94e−52), CD4+T cells (Cor=0.299, p=1.19e−10), macrophages (Cor = 0.247, p = 9.59e−8), neutrophils (Cor = 0.612, p = 7.03e−48), and dendritic cells (Cor = 0.648, p = 1.56e−54; Fig. 10a). Similar results were obtained for CXCL10, CXCL13, CCL4 and CCL5 (Fig. 10b-e). While, the correlation between NMU/GAL and immune infiltration is not obvious (Fig. 10f-g). Additionally, we also found significant correlations of these seven genes with 28 types of TILs across human heterogeneous cancers (Additional file 2: Fig. S1). CXCL9 significantly positive correlated with abundance of 28 types TILs such as activated CD8 T cells (Act_CD8 T cells; rho = 0.801, p < 2.2e-16), regulatory T cells (Treg; rho = 0.693, p < 2.2e-16), macrophage (rho = 0.655, p < 2.2e-16), natural killer T cells (NK T cells; rho = 0.74, p < 2.2e-16), myeloid derived suppressor cells (MDSC; rho = 0.744, p < 2.2e-16), immature B cell (Imm_B; rho = 0.78, p < 2.2e-16). Similar results were found in CXCL10, CXCL13, CCL4 and CCL5 (Additional file 2: Fig. S1).
Additionally, in the survival model of TIMER, results showed that high levels of B cells (Log-rank P=0), CD8+T cells (Log-rank P=0), neutrophils (Log-rank P=0), dendritic cells (Log-rank P=0), CXCL9 (Log-rank P=0), CXCL10 (Log-rank P=0), CXCL13 (Log-rank P=0), CCL4 (Log-rank P=0), CCL5 (Log-rank P=0), were significantly related to longer survival time in SKCM patients (Additional file 3: Fig. S2).
Drug–gene interaction prediction
We obtained 42 drug-gene interaction pairs in DGIdb, including genes (IL6, CXCL2, CXCL10, CCL4, CCL5, PTGER3, GAL, NPY1R) and 69 drugs, as shown in Fig. 11. This result may help develop new targets for SKCM therapy.