3.1. Patient characteristics
In February 2020, the gene expression profiles and relevant clinicopathological variables of 471 patients with SKCM were extracted from the TCGA database. Among these patients, 293 were male (62.21%) and 178 were female (37.79%). The median age of the included cases was 58 years (range, 15–90 years). In this study, the numbers of patients with tumor grades G1, G2, G3, and G4 were 5 (1.1%), 18 (3.8%), 77 (16.3%), and 167 (35.5%), respectively. In total, 78 cases were diagnosed with stage I SKCM (16.6%), 141 were diagnosed with stage II SKCM (29.9%), 175 were diagnosed with stage III SKCM (37.2), and 23 were diagnosed with stage IV SKCM (4.9%). The most common tumors were in T4 group (33.1%, n = 156); T1, T2 and T3 groups were 8.9% (n = 42), 16.3% (n = 77) and 20.0% (n = 94), respectively. In addition, 184 patients had regional lymph node metastasis while 24 had distant metastasis.
3.2. Correlation between immune/stromal scores and clinicopathologic variables of skin cutaneous melanoma
In this study, complete gene expression profiles and clinicopathological information from the TCGA database were used to analyze the data of all cases. According to the data of the ESTIMATE algorithm, stromal scores ranged from −1806.74 to 1892.23, immune scores ranged from −1481.02 to 3769.31, and ESTIMATE scores ranged from −2632.50 to 5077.47. The mean immune (P=0.009), stromal (P=0.012), and ESTIMATE (P=0.006) scores of patients with SKCM in the older-age group were consistently lower than those of patients in the lower-age group (Fig. 1A). Moreover, the mean immune (P=1.78e−04), stromal (P=3.69e−03), and ESTIMATE (P=3.95e−04) scores differed among stages I, II, III, and IV(Fig. 1B); here, the overall trend showing an initial decrease followed by a gradual increase. The mean immune (P=5.14e−04), stromal (P=0.017), and ESTIMATE (P=1.83e−03) scores also varied according to tumor stage(Fig. 1C). These results suggest that stromal, immune, and assessment scores are correlated with certain clinicopathologic variables.
Kaplan–Meier survival curves were drawn to assess the potential associations between overall survival and immune, stromal, and assessment scores. First, the 471 patients with SKCM were divided into high and low groups according to their scores. Immune score results showed that, compared with the high-score group, the overall survival time of patients with SKCM in the low-score group was longer (Fig. 2A, P=0.001). ESTIMATE score results showed that the overall survival time of patients with SKCM in the low-score group was longer than that in the high-score group (Fig. 2, P=0.001). Finally, stromal scores (Fig. 2B, P=0.069) indicated that the overall survival time of patients in the low-score group was longer than that in the high-score group, but the difference was not statistically significant.
Next, we analyzed Affymetrix chip data from our 471 patients to reveal differences in gene expression profiles based on immune and stromal scores. A heat map drawn on the basis of immune scores is shown in Fig. 3A, which showed 1049 and 94 up-regulated genes in the high- and low-score groups, respectively (|fold change|>2, P<0.05). Fig. 3B shows the heatmaps drawn on the basis of stromal scores; here, 1369 and 30 genes were up-regulated in the high- and low-score groups, respectively (|fold change|>2, P<0.05). Differentially expressed genes (DEGs) in the immune and stromal score groups were identified by Venn map. Our results showed that 909 genes were up-regulated in the high-score group (Fig. 4A) while 5 genes were down-regulated in the low-score group (Fig. 4B).
3.3. Functional enrichment analysis of differentially expressed genes
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) were used for functional enrichment analysis to determine the functions of DEGs. GO analysis was employed to determine the top 10 biological processes (BP), top 10 molecular functions (MF), and top 10 cellular components (CC) (Fig. 5A). The top 10 enriched BPs included leukocyte cell–cell adhesion, T-cell activation, T-cell aggregation, lymphocyte aggregation, leukocyte aggregation, adaptive immune response, regulation of leukocyte activation, and regulation of lymphocyte activation. These BPs are involved in the tumor immune response. The top 10 enriched MFs included cytokine activity, cytokine receptor activity, chemokine activity, cytokine receptor binding, chemokine receptor binding, CCR chemokine receptor binding, antigen binding, carbohydrate binding, G-protein coupled chemoattractant receptor activity, and chemokine receptor activity. These MFs are also associated with the tumor immune response. The top 10 CCs included side of membrane, external side of plasma, MHC class II protein complex, T-cell receptor complex, MHC protein complex, receptor complex, immunological synapse and plasma membrane receptor Integral component of luminal side of endoplasmic reticulum. KEGG analysis yielded similar results (Fig. 5A). KEGG pathway analysis also indicated that the cytokine–cytokine receptor interaction and chemokine signaling pathways are involved in the tumor immune response and SKCM development.
The STRING tool was used to explore the protein–protein interaction (PPI) network and elucidate interactions between DEGs. The functional network consisted of 14 modules, 339 nodes and 848 edges (Fig. 6A). The first 20 significant genes were CCR2, CD8B, CXCL11, CXCL13, CXCR3,CD72, CCL25, CD3G, C1R, C1S, CXCL9, CCL21, CCR7, CXCR5, CCR8, CXCL16, CXCL10, C1QC, HLA-DRB5 and CCL19 (Fig. 6B).
3.4. Role of individual DEGs in overall survival in skin cutaneous melanoma
In this study, Kaplan–Meier survival curves were drawn to determine the potential function of 909 up-regulated and 5 down-regulated genes in the overall survival of SKCM patients. The log-rank test showed that, among 1014 DEGs, AQP9, CCR8, PATL2, RGL4, OR9G1, CXCR5, IPCEF1, FCRL3, RCSD1, CP, NEXN, MAFB, SLC40A1, and TNFRSF13B have significant predictive effects on the overall survival of SKCM ( Fig. 7, all P<0.05).
3.5. Identification of the kinase target of hub genes
Among the top 20 important genes (i.e., CD72, CXCR5, CCR8, CXCL13, CXCR3, CXCL16, CXCL10, CCR2, CD8B, CXCL11, CCL25, CD3G, C1R, C1S, CXCL9, CCL21, CCR7, C1QC, HLA-DRB5, and CCL19 ) found in the PPI network, CCR8 and CXCR5 were significantly correlated with the overall survival of SKCM patients. Therefore, these DEGs were selected as hub genes. We explored the kinase targets of these genes by using the LinkedOmics database and found that PRKCA, ATM, PLK1, CHEK1, and ATR are the five most significant kinase target networks for CXCR5. Moreover, kinases LYN, FYN, AURKB, PAK1, and LCK were the five most prominent target networks for CCR8 (Table 1). We also studied the co-expression genes of CXCR5(Fig. 8A) and CCR8 (Fig. 8B).
3.6. Immune landscape of the tumor microenvironment in melanoma
We analyzed the correlation between hub genes and immune cell infiltration. Our results revealed a positive correlation between CCR8 expression and CD8+ T cells (Cor=0.592, P=8.51e−43), CD4+ T cells (Cor=0.389, P=1.34e−17), B-cell infiltration (Cor=0.327, P= 1.17e−12), macrophages (Cor=0.412, P=5.83e−20), neutrophils (Cor=0.662, P=2.27e−58), and dendritic cells (Cor=0.625, P=1.33e−49) (Fig. 9A). Similar results were also obtained for CXCR5. Specifically, CXCR5 expression was positively correlated with CD8+ T cells (Cor=0.413, P=1.63e−19), CD4+ T cells (Cor=0.49, P=2.48e−28), B-cell infiltration (Cor=0.419, P=1.75e−20), macrophages (Cor=0.223, P=1.62e−06), neutrophils (Cor=0.398, P=1.21e−18), and dendritic cells (Cor=0.488, P=4.02e−28) (Fig. 9B).