Identification of DEGs
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 (Figure 2a) (Additional file 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, displayed in bubble charts (p<0.05, Figure 3). GO analysis indicated that changes in biologic processes of DEGs were significantly enriched in the inflammatory response, immune response, regulation of transcription. In addition to cellular component, the overlapping 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. At last, 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 demonstrated a correlation between DEGs and immune response and inflammatory response.
PPI network construction and hub genes analysis
The PPI network was constructed via Cytoscape, including 247 nodes and 792 edges (Figure 2b). Then, the most important PPI network module was obtained using MCODE, consisted of 12 nodes and 64 edges (Figure 2c). After selection, CCL5, PTGER3, IL6, CXCL13, CCL27, CCL4, CXCL9, NMU, CXCL2, GAL, NPY1R, CXCL10 were considered as hub genes of the network. The pathway analysis of hub genes via ClueGO (Figure 4) revealed that the hallmark pathways include interleukin-10 signaling, peptide ligand-binding receptors and chemokine receptors bind chemokine signaling.
Prognostic and diagnostic value assessment of hub genes
We found that CCL4, CCL5, CXCL9, CXCL10, CXCL13 were significantly positively associated with overall survival (OS) and NMU, GAL were negatively related to OS (P<0.01) (Figure5). In addition, CCL4(AUC=0.872), CCL5(AUC=0.775), CXCL9(AUC=0.834), CXCL10 (AUC=0.786), CXCL13 (AUC=0.807), NMU(AUC=0.829) except GAL(AUC=0.536) were of diagnostic value with p value<0.001 (Figure 6).
Validation of prognostic genes expression and clinical characteristics
In the present study, 472 melanoma patients from the TCGA cohort were used to analyze correlation between prognostic gene expression and pathological stages, Breslow depth, which are useful prognostic indicators for melanoma. As shown in figure 7, CXCL9, 10, 13 and CCL4, 5 were significantly overexpressed in stage Ⅰ and decreased in subsequent stages with P value<0.001. What’s more, the expressions of CXCL9, 10 and CCL4 (p<0.05) were negatively correlated with Breslow depth from (0-1.5mm) to (>3mm) while NMU, GAL (P<0.05) were positively correlated with it (Figure 8). Moreover, CXCL9, 10, 13 and CCL4, 5, NMU were high expressed in tumor tissue compared to normal (p<0.01, Figure 9). IHC images suggested that CCL4, 5, CXCL9,10, 13 and NMU were detected in melanoma tissues while not detected in normal tissues (Figure 10). All the results suggested these hub genes are of important prognostic value.
Immune infiltration analysis
After determining the prognostic value, we performed correlation analysis between prognostic genes and immune infiltration level for melanoma. 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; Figure 11a). Similar results were obtained for CXCL10, CXCL13, CCL4 and CCL5 (Figure 11b-e). While GAL was found negatively correlated with infiltration of CD8+T cells (Cor= −0.171, p=3.22e−4), CD4+T cells (Cor= −0.108, p=2.30e−2), neutrophils (Cor = −0.185, p = 7.94e−5), and dendritic cells (Cor = −0.135, p = 4.42e−3; Figure 11g). Additionally, we also found significant correlations of these seven genes with 28 types of TILs across human heterogeneous cancers (Figure 12). 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 (lmm_B; rho = 0.78, p < 2.2e-16) in Figure 12a. Similar results were found in CXCL10, CXCL13, CCL4 and CCL5 (Figure 12b-e). While NMU was negatively correlated with abundance of TILs (Act_CD8 T cells; rho = −0.159, p < 0.001, Treg; rho = −0.138, p < 0.01, macrophage; rho = −0.108, p < 0.05, NK T cells; rho = −0.176, p < 0.001, MDSC; rho = −0.094, p <0.05, lmm_B; rho = −0.207, p <0.001) and similar results were found in GAL as shown in Figure 12f-g.
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 Figure 13. This result may help to develop new targets for melanoma therapy.