Skip to main content

Gene expression and immune infiltration in melanoma patients with different mutation burden



Immunotherapy is a vital component in cancer treatment. However, due to the complex genetic bases of cancer, a clear prediction index for efficacy has not been established. Tumor mutation burden (TMB) is one of the essential factors that affect immunotherapeutic efficacies, but it has not been determined whether the mutation is associated with the survival of Skin Cutaneous Melanoma (SKCM) patients. This study aimed at evaluating the correlation between TMB and immune infiltration.


Somatic mutation profiles (n = 467), transcriptome data (n = 471), and their clinical information (n = 447) of all SKCM samples were downloaded from The Cancer Genome Atlas (TCGA) database. For each sample, TMB was calculated as the number of variants per megabase. Based on K-M survival analysis, they were allocated into the high-TMB and low-TMB groups (the optimal cutoff was determined by the ‘surv_cutpoint’ algorithm of survival R package). Then, Gene ontology (GO) and Gene Set Enrichment Analyses (GSEA) were performed, with immune-associated biological pathways found to be significantly enriched in the low-TMB group. Therefore, immune genes that were differentially expressed between the two groups were evaluated in Cox regression to determine their prognostic values, and a four-gene TMB immune prognostic model (TMB-IP) was constructed.


Elevated TMB levels were associated with better survival outcomes in SKCM patients. Based on the cutoff value in OS analysis, they were divided into high-TMB and low-TMB groups. GSEA revealed that the low-TMB group was associated with immunity while intersection analysis revealed that there were 38 differentially expressed immune-related genes between the two groups. Four TMB-associated immune genes were used to construct a TMB-IP model. The AUC of the ROC curve of this model reached a maximum of 0.75 (95%CI, 0.66–0.85) for OS outcomes. Validation in each clinical subgroup confirmed the efficacy of the model to distinguish between high and low TMB-IP score patients.


In SKCM patients, low TMB was associated with worse survival outcomes and enriched immune-associated pathways. The four TMB-associated immune genes model can effectively distinguish between high and low-risk patients.

Peer Review reports


Skin Cutaneous Melanoma is a fatal malignant tumor of the skin. Currently, the specific etiology of primary cutaneous melanoma is postulated to be long-term ultraviolet overexposure [1, 2]. Ultraviolet short waves induce DNA damage in skin cells, which leads to inflammation, immunosuppression, and melanoma [3]. Even though SKCM accounts for less than 5% of skin cancers, it has a high mortality rate that is correlated with its high malignancy and invasiveness [4]. Over the past 30 years, the global prevalence of SKCM has been on the rise, and it has been reported that in 2018, global incidences were about 280 thousand new morbidites and over 60 thousand mortalities [5]. In 2019, there were more than 96 thousand new cases of SKCM in the USA alone [6]. In China, the annual number of new SKCM cases is more than 8 thousand [7]. Even though SKCM incidences are increasing, advances in medical technology, especially in targeted therapy, immunotherapy, and chemotherapy to a certain extent, have led to positive survival outcomes for patients and reduced mortality rates. For example, before 2011, the median survival time for metastatic melanoma patients was 9 months, but has now been reported to exceed 2 years [8]. These improvements are mainly attributed to small molecule inhibitors (e.g., BRAF inhibitors, MEK inhibitors) [9, 10] and immune checkpoint blockade (ICB) [11,12,13]. Identification of CTLA-4 and PD-L1/PD-1 antibodies has enhanced advances in tumor immunotherapy [14, 15].

The tumor microenvironment is mainly composed of tumor cells, fibroblasts, immune cells and the extracellular matrix, which significantly affect treatment and survival outcomes. Immunotherapy is closely correlated with immune cells. Elucidating on immune infiltration in the tumor microenvironment is key to improving response rates, and could inform the development of new immunotherapeutic strategies. Melanoma is a tumor with strong immunogenicity. Pathological studies have reported that there is a high number of immune cell infiltration in melanoma tissues and, immunotherapy can directly inhibit tumor progression or even cure tumors [16, 17].

Tumor mutation burden refers to the frequency of mutations in the coding regions of somatic cells (variants per megabase) [18,19,20]. An increase in TMB is correlated with an increase in tumor antigenicity, which is the premise of the effectiveness of the PD-1/PD-L1 antibody. In several clinical trials, TMB has been shown to be a good predictor of immunotherapeutic efficacy [21,22,23,24,25]. The CheckMate-026 clinical trial, a retrospective analysis, reported that: among NSCLC patients administered with immunotherapy, the remission rate and progression-free survival outcomes of the high-TMB group were significantly better than those of the low-TMB group. Patients exhibiting elevated PD-L1 protein expression with high TMB had the most significant benefits, however, patients with elevated PD-L1 protein expression levels but with low TMB had no significant benefits. CM012 post-test analysis revealed corresponding results [21].

Various public databases, such as TCGA, GEO, cBioPortal, and ICGC among others, are available for researchers. However, few studies have evaluated TMB associated immune infiltration in SKCM. Therefore, this study aimed at evaluating the prognostic value of TMB, and to elucidate on how it is associated with immune infiltration.


Data acquisition

Somatic mutation profiles for 467 SKCM patients were retrieved from the TCGA database ( The “Maftools” R package [26] was used to represent the mutation situation. Moreover, we downloaded the level 3 transcriptome data for all available SKCM samples (tumor samples, n = 471). Corresponding clinical information including sex, age, TNM stages, pathological stage, as well as survival outcomes were also obtained (n = 447, Table 1, Supplementary Table S1). Samples with follow-up time of less than 60 days were deleted after which the remaining samples were merged with TMB for survival analysis (Supplementary Table S2). The workflow of this study was illustrated in Supplementary Fig. S1. These data were retrieved from free public databases, and as such, ethical approval was waived.

Table 1 Clinical baseline of 447 SKCM patients

TMB calculation and Kaplan-Meier analysis

TMB refers to the total number of substitutions, insertions, deletions, and mutant genes per megabase in the coding region (exon) of the gene assessed in the tumor tissue. In this study, we determined TMB by dividing the number of variants by the length of exons (38 million) for each sample. TMB and the corresponding survival time of the same sample were merged. Subsequently, Kaplan-Meier (KM) analysis was performed to compare survival outcomes in low- versus high-TMB groups and determined the p of the log-rank test. Moreover, Wilcoxon rank test was performed to assess the difference between two groups of different clinical characteristics while the Kruskal-Wallis test was performed to compare differences among multiple groups.

Differentially expressed genes and functional pathway analysis

RNA-seq data of SKCM patients were divided into low- and high-TMB groups. The “limma” R package [27] was used to identify DEGs between the two groups (Fold Change [FC] = 1 and False Discovery Rate [FDR] < 0.05)., and the Heatmap plot was drawn using “pheatmap” R package ( to indicate the difference. Next, the Entrez ID for every DEG was obtained using the “org.” package, after which GO analysis was performed using “clusterProfiler” [28], “enrichplot” and “ggplot2” packages. In addition, GSEA [29] was performed using the TMB level as the phenotype. The “c2.cp.kegg.v7.1.symbols.gmt” was retrieved from the MSigDB database [30], and was set as a baseline gene set. Pathways with FDR < 0.25 were considered enriched. Due to differences in immune pathway enrichment in high versus low TMB groups, we used the “estimate” R package ( to calculate immune cell scores for the transcriptome data. From the Immunoscore survival analysis, it was shown that low immune infiltration is associated with poor survival outcomes. Therefore, interactions between DEGs and immune-related genes (ImmPort Private Database) were evaluated. Venn analysis showed that there were 38 intersection genes for 1812 immune genes [31] and 504 DEGs.

Construction of a TMB-immune prognostic model for hub immune genes

Univariate-lasso-multivariate Cox regression analysis was performed on the 38 intersection genes, from which four hub immune genes were obtained and used to construct the TMB-IP model. Then, we calculated the TMB-IP score for all patients by the coefficients of each gene and divided the SKCM patients into high and low TMB-IP group. ROC curve and K-M analysis were used to assess the predictive value of the TMB-IP score in SKCM. The prognostic efficiency of this model was also tested in each clinical subgroup through K-M survival analysis.

TIMER database and CIBERSORT algorithm

Mutation types of the hub immune genes with different immune infiltrates in SKCM were assessed according to the “SCNA” module of Tumor Immune Estimation Resource (TIMER) database [32]. Furthermore, “CIBERSORT” algorithm [33] was used to estimate the immune infiltration degree of 22 types of cells in a mixed population of cells based on certain features of gene expression in 22 leukocyte subtypes-LM22. The “pheatmap” package revealed immune cell distributions in the two groups. Then, we used the Wilcoxon rank test to assess disparities in the amounts of immune infiltrates in the low-versus high-TMB groups. The “vioplot” R package was used to determine the p-values.

Determination of prognostic value of immune cells in the TIMER database

Data from the TIMER database were used to perform multivariate Cox analysis of the cells that were involved in immune infiltration, as well as to calculate the hazard ratio (HR; 95%CI). Survival outcomes were evaluated by survival analysis.

Statistical analysis

R software (Version 3.5.2) was used to conduct all the data analysis in the present study. The optimal cutoff for samples classified into better and worse survival groups was determined by the ‘surv_cutpoint’ algorithm of survival R package.


Mutation profiles in SKCM

Somatic mutation profiles for 467 SKCM patients were retrieved from the TCGA database. According to the mutation data, “Maftools” R package was used to examine the findings. The waterfall plot was used to present the mutation data for every gene in every sample (Fig. 1a). Further, mutations were grouped based on various categories. In the grouping, missense mutations were the most common (Fig. 1b), while deletion/insertion mutaions were common than single nucleotide polymorphisms (SNP) (Fig. 1c). Regarding single nucleotide variants in SKCM, C > T occurred more often (Fig. 1d). In addition, we constructed an SKCM box plot (different colors denote different mutations) to reveal mutation types based on the number of changed bases in every sample (Fig. 1e, f). Figure 1g shows the top ten mutated genes in SKCM, which were TTN (72%), MUC16 (67%), BRAF (51%), DNAH5 (49%), PCLO (44%), LRP1B (38%), ADGRV1 (35%), RP1 (33%), ANK3 (32%) and DNAH7 (32%). Supplementary Fig. S2 shows coincidences, as well as exclusive relationships among the mutated genes (green denotes co-occurrence, whereas red denotes relationships that are mutually exclusive). It is shown that genes with higher mutations in SKCM appear at the same time, and that there are no obvious repulsive genes. The genecloud plot was established to show the frequency of mutations in other genes (Supplementary Fig. S2).

Fig. 1

Mutation profile landscape in SKCM samples. a Waterfall plot showing the mutation details of every gene in every sample; b,c,d Various types of mutation classified based on different groups; e, f burden of tumor mutation in particular samples; g Top ten mutated genes. Various types of mutations are represented by different color annotations at the right-bottom

TMB is correlated with survival outcomes, pathological stage and tumor grade

Based on the above mutation data profile, we computed the number of mutation events/million bases as TMB for the 467 samples. The optimal cutoff for samples classified into better and worse survival groups was established by the ‘surv_cutpoint’ algorithm of survival R package, consistent with the recognition that higher TMB enhances immune recognition and leads to better disease outcomes. High-TMB group patients exhibited better OS outcomes (p < 0.001, HR = 0.54, 95%CI = 0.39–0.75; Supplementary Table S2, Fig. 2a), better DSS (p < 0.001, HR = 0.49, 95%CI = 0.35–0.7; Supplementary Table S2, Fig. 2b) and better PFI (p = 0.003, HR = 0.68, 95%CI = 0.52–0.9; Supplementary Table S2, Fig. 2c) than the low-TMB group. Moreover, elevated TMB levels were correlated with advanced age (p = 0.0045; Fig. 2d), male (p = 8.5e− 05; Fig. 2e), lower tumor pathological stages (stage I + II vs. stage III + IV, p = 0.022; Fig. 2f) and lower AJCC-N stages (N0 vs. N3, p = 0.002; N2 vs. N3, p = 0.034; Fig. 2h). However, there were no significant differences between the associations of TMB and AJCC-T (Fig. 2g) and AJCC-M (Fig. 2i).

Fig. 2

TMB prognosis and its relationship with clinical features. Higher TMB levels correlated with better OS (a), better DSS (b), and better PFI (c); d-i Wilcoxon and Kruskal-Wallis test in groups of different clinical characteristics

Differential expression analysis between two groups

Among the 471 RNA-seq samples, a total of 467 samples were selected from the same source as mutation data. Based on the TMB cutoff value of the OS analysis (cut-off = 6.37, Fig. 2a), we divided the RNA-seq data for 467 tumor samples into low (n = 148) and high-TMB (n = 319) groups (Fig. 3a). Using the “limma” R package, a total of 504 DEGs were identified between these two transcriptomic data sets with FC = 1 and FDR < 0.05. GO enrichment analysis revealed that the DEGs were mainly involved in epidermis development, skin development, and other biological processes (Fig. 3b). Further, we analyzed the GSEA data of the top TMB-associated items. Figure 3c shows that the low-TMB group was associated with immunity, and was mainly concentrated in autoimmune thyroid disease, B and T cell receptor signaling pathways, as well as intestinal immune network for the production of IGA, while the high-TMB group lacked this function (Supplementary Table S34). Given that TMB was associated with immune signature/pathways in SKCM, various analyses were performed on the tumor microenvironment. Using the transcriptome data, the “estimate” R package was used to calculate immunescores for all samples (Non-tumor samples were deleted, OS time over 60 days were merged, n = 442, Supplementary Table S5). The score represents the infiltration degree of immune cells in each sample. The higher the score, the more immune cells infiltrated the sample. Immunoscore survival analysis revealed that lower immune infiltration is associated with poor prognosis (p < 0.001, HR = 0.47, 95%CI = 0.36–0.63; Fig. 3d). These findings imply that the lower TMB led to immune-associated pathway enrichment while the lower immune infiltrate group was correlated with poor prognosis. Therefore, we evaluated the predictive accuracy of the TMB-associated immune genes. Venn analysis showed that, from 1812 immune genes and 504 DEGs, there were 38 intersecting genes (Supplementary Fig. S3).

Fig. 3

Differentially expression analysis and Construction of TMB-IP. a Top 50 DEGs were displayed in heatmap plot; b GO results of 504 DEGs; c GSEA revealed the top TMB- associated pathways revealed that low-TMB was correlated with immune-related pathways; d Survival analysis based on the immunescore; e-g Univariate-lasso-multivariate Cox regression of 38 immune-related DEGs; h Survival analysis based on the TMB-IP score; i Multi-year ROC curves based on the TMB-IP score; j-m Survival analysis of four immune-related hub genes based on their RNA expression levels; n-s Survival analysis of six clinical subgroups based on high and low TMB-IP grouping

Construction and assessment of TMB-IP in SKCM

Prognostic values of the 38 immune-related DEGs were further evaluated. Data for the 442 SKCM patients were used in univariate-lasso-multivariate Cox regression analysis. A total of 11 genes obtained from the univariate Cox analysis (p < 0.001; Fig. 3e) while 7 genes were obtained from lasso regression analysis (PGLYRP3, GAL, ADCYAP1R1, SLPI, CNTFR, LIF, and CRABP2; Fig. 3f). Finally, PGLYRP3, GAL, ADCYAP1R1 and LIF (Fig. 3g) were selected to construct the Tumor Mutation Burden Immune Prognostic model (TMB-IP) using the following formula: TMB-IP score = (PGLYRP3*0.0988045135949462 + GAL*0.0483710471635067 + ADCYAP1R1*0.123387518699318 − LIF*0.121787959868336). We calculated the TMB-IP score for each SKCM patient (Supplementary Table S6). The K-M plot revealed that high TMB-IP score patients exhibited worse survival outcomes (Fig. 3h). Using a cutoff of 0.97080003, patients were divided into high (n = 231) and low (n = 211) TMB-IP groups. We drew the ROC curves to evaluate the predictive accuracyof the model. The OS outcomes were: AUC = 0.69 (95%CI = 0.60–0.77) for 1-year, 0.68 (95%CI = 0.63–0.73) for 5-years, 0.69 (95%CI = 0.63–0.75) for 10-years, 0.73 (95%CI = 0.64–0.82) for 15-years, 0.75 (95%CI = 0.66–0.85) for 17-years and 0.72 (95%CI = 0.58–0.85) for 20-years (Fig. 3i). Batch survival analysis indicated that these four prognostic hub immune genes were closely correlated with survival outcomes. Suppressed expression levels of ADCYAP1R1 (p < 0.001, HR = 1.75, 95%CI = 1.32–2.3; Fig. 3j), GAL (p < 0.001, HR = 1.94, 95%CI = 1.48–2.55; Fig. 3k) and PGLYRP3 (p < 0.001, HR = 1.67, 95%CI = 1.24–2.25; Fig. 3m) were positively correlated with better prognosis, while suppressed expression levels of LIF (p < 0.001, HR = 0.55, 95%CI = 0.42–0.73; Fig. 3l) were negatively correlated with better survival outcomes. Importantly, we further assessed the prognostic efficiency in different clinical subgroups. It was found that high and low TMB-IP score classification criteria can perfectly distinguish the prognostic outcomes for: female (p < 0.001, HR = 0.32, 95%CI = 0.19–0.54; Fig. 3n), male (p < 0.001, HR = 0.5, 95%CI = 0.35–0.7; Fig. 3o), stage I + II (p < 0.001, HR = 0.37, 95%CI = 0.24–0.56; Fig. 3p), stage III + IV (p < 0.001, HR = 0.43, 95%CI = 0.28–0.68; Fig. 3q), ≥60Y (p < 0.001, HR = 0.51, 95%CI = 0.34–0.77; Fig. 3r) and < 60Y (p < 0.001, HR = 0.35, 95%CI = 0.24–0.52; Fig. 3s) patient groups.

Relationships between mutants and immune infiltrates

The link between the mutants of the 4 hub genes and immune infiltrates in the SKCM microenvironment was assesssed. Based on the findings, immune infiltrate inhibitions, such as CD8+ T cells, dendritic cells, CD4+ T cells, neutrophils, B-cells, and macrophages depended on the type of mutation exhibited by the genes, relative to the levels of immune infiltration in the wild type samples (Fig. 4).

Fig. 4

Associations mutants with immune cells infiltration of four hub genes. a-d Mutants (Arm-level Deletion, Deep Deletion, High Amplication, Arm-level Gain) of four TMB-associated genes exhibited low level of immune cell infiltration compared with Diploid/Normal

Variations in the abundance of immune cell infiltration in the low- versus high-TMB groups

We have shown that the DEGs negatively impacted on immune pathways and that mutants of the four hub genes were inversely correlated with immune infiltrates. We further compared the various immune fraction profiles in the high- versus low-TMB groups. We used the Voom algorithm in “limma” R package to normalize the transcriptome data for SKCM samples. Samples were filtered at p > 0.05 using the “CIBERSORT” R package, and 201 samples (125 low-TMB samples and 76 high-TMB samples) were identified, which were then used to analyze the immune cells (Supplementary Table S7). A box plot was constructed to show specific fractions of 22 immune cells based on “CIBERSORT” algorithm in every SKCM sample (Fig. 5a). These 22 types of immune cells are subclasses of the 6 immune cells in Fig. 4. Wilcoxon rank test revealed that infiltration levels of CD8+ T cells, CD4+ memory activated T cells, follicular helper T cells, monocytes and macrophage M1 cells were reduced in the high-TMB group, relative to those in the low-TMB group (Fig. 5b). Correlation and distribution analysis of the infiltration degree of 22 immune cells in the two groups are shown in Supplementary Fig. S4.

Fig. 5

Twenty-two immune fractions in low vs. high-TMB groups. a Barplot of the 22 immune fractions denoted by different colors: b Wilcoxon test of 22 immune fractions in high- vs. low-TMB group

Association between immune cell infiltration and survival

The Cox regression model was used to determine the association between immune cells and prognosis in SKCM samples (Table 2). In the model, Survival (SKCM) = B cell + CD8 T cell + CD4 T cell + Macrophage + Neutrophil cell + Dendritic cell + ADCYAP1R1 + LIF + GAL + PGLYRP3. Elevated B cell, CD4+ T cell, Macrophage cell infiltrates eas well as elevated expression levels of ADCYAP1R1, GAL, PGLYRP3 were found to be risk factors for SKCM (HR > 1). Additionally, we performed the K-M analysis (Supplementary Fig. S5), where elevated infiltration levels of B cells (P = 0), CD8+ T cells (P = 0), Neutrophils (P = 0) and Dendritic cells (P = 0) were positively correlated with better SKCM prognosis, consistent with the result of the survival analysis of the overall infiltration of immune cells in the microenvironment (Fig. 3).

Table 2 Multivariate Cox regression analysis of six immune infiltration cells and four immune-related genes


We analyzed the SKCM-cohort in TCGA. It was found that higher TMB is associated with better survival outcomes and that transcriptomic data for low-TMB patients was enriched in immune-related pathways. Prognostic model construction confirmed that the four immune-related genes can effectively distinguish patients with different prognostic outcomes.

Melanoma, a highly proliferating tumor, originates from melanocytes in the neural crest. It occurs in the skin, intestinal tract, and genital mucosa. Combinations of immunotherapy and targeted therapy has completely changed the treatment of SKCM. In 2013, Science ranked tumor immunotherapy at the top of ten scientific breakthroughs [34]. Immunotherapeutic agents targeting CD8+ T-cell surface receptors (CTLA-4 and PD1) have greatly improved patient prognosis [8]. Snyder et al. found that TMB is correlated with clinical benefits of immunotherapy, that is, the greater the number of somatic mutations, the more likely the tumor is to respond to ICB [35,36,37].

Through integration and unified processing of high-throughput sequencing data for multiple samples from the TCGA database, we established significant mutation genes that may be involved in immune responses of SKCM. It is important to establish molecular markers of immunotherapy to inform individualized treatment of melanoma patients using immunotherapy.

Survival analysis revealed that the high-TMB group patients exhibited better survival outcomes than those in the low-TMB group. Median survival outcomes for the high-TMB group and low-TMB group were 1438 days and 698 days, respectively. Melanoma is a tumor type with a high number of immune cells. Some mutations produce new antigens that are recognized by T cells, which may be one of the reasons for the better survival outcomes of the high-TMB group. Besides, higher TMB correlated with advanced age, male, lower tumor pathological stages, and lower AJCC-N stages. These findings are consistent with those of other studies. Rizvi et al. reported that NSCLC patients with elevated non-synonymous mutation burdens can achieve long-term remission by immunotherapy, which is of paramount importance for PFS [38]. The CheckMate-227 trial, comparing nivolumab+ipilimumab, nivolumab, and chemotherapy, revealed that patients with high TMBs (> 10 mutations/Mb) have the best OS outcomes [39], supporting the validity of TMB-based treatment stratification.

Differential expression analysis and GSEA showed that the low-TMB group was correlated with autoimmune thyroid disease, B-and T-cell receptor signaling pathways, and intestinal immune networks for IGA production. This finding is unexpected and surprising because there is no enrichment of any immune-related pathway in the high-TMB group. GSEA is a method for enriching biological pathways based on hundreds of gene sets. The close relationships between low-TMB and multiple immune-related pathways have been reported, but the specific mechanisms should be further evaluated. Moreover, survival analysis based on immune score of microenvironments indicated that higher levels of immune infiltration are associated with better survival outcomes. Infiltration levels of a variety of immune cells, mainly T cells, can react with tumor or non-tumor antigens, which affects survival outcomes for patients. We also analyzed the effect of each type of immune cell infiltration on survival.

Finally, a four-immune gene prognostic model was developed based on 38 intersecting genes from 1812 immune genes and 504 DEGs, which are greatly useful for survival prediction. Four hub genes were either positively correlated (LIF) or negatively correlated (PGLYRP3, GAL, and DCYAP1R1) with prognosis. The median survival outcomes for high TMB-IP score and low TMB-IP score groups were 804 days and 1490 days, respectively. TMB-IP score is a risk score of the prognostic model. High TMB-IP score means high risk, and it has a relatively shorter median survival outcomes than the low TMB-IP score (low risk) group. Similar findings were obtained in the survival analysis of the six subgroups.

Furthermore, mutants of these hub immune genes were associated with immune infiltrates, such as CD8+ T cells, CD4+ T cells, dendritic cells, macrophages, and B-cells in the SKCM microenvironment, implying that these immune infiltrate cells might be suppressed by mutations. To some extent, this is consistent with the finding of Rosenthal: the absence of new antigens in high immune infiltrated lung cancer suggests that these tumors may evade immune attacks by inhibiting the expression of new antigens [40]. Differences in immune infiltration may affect tumor immune editing and the appearance of neoantigens in tumors. During tumor evolution, tumors with low immune infiltrations exhibit reduced tumor antigen editing, indicating that the historical immune editing or the original clone neoantigen copy is lost. It has been found that immune cell activation in the tumor microenvironment is not the same in different tumor types, even in the same mutation feature, therefore, a specific analysis is needed [41, 42]. We have shown that in SKCM patients, infiltration levels of CD4+ memory activated T cells, CD8+ T cells, follicular helper T cells, Monocytes, and Macrophage M1 cells in the high-TMB group were significantly suppressed. Even though there were no significant differences in other categories of immune cells, the degree of infiltration was not an absolute trend as described above. There were immune cells with higher infiltration levels in the high-TMB group, and specific reasons for this outcome are worthy of further investigations.

In the tumor microenvironment, different immune cells confer different tumor responses. T-helper-1 cells, natural killer cells, M1 phenotype macrophages, and DC1 phenotype dendritic cells are involved in the suppression of tumorigenesis and development, while T-helper-2 cells, M2 macrophages, DC2 dendritic cells, and regulatory T cells (Treg) suppress immune responses [43,44,45]. Analysis of the unique properties of immune cells in the tumor microenvironment may inform the design of cancer immunotherapy targets. The correlation between lymphocyte infiltration in the tumor microenvironment and immunotherapeutic benefits have been confirmed [46,47,48,49]. Higher infiltration levels of B cells, CD8+ T cells, Neutrophils, and Dendritic cells are positively correlated with better SKCM prognosis. Changes in tumor immune microenvironment, tumor gene mutation, and gene regulation may affect tumor evolution and survival outcomes.

However, there are some imperfections. First, we did not carry out experiments (vivo/vitro) to verify the relationship between four immune genes and immune infiltration. Second, the sample size was not adequate to confirm the potential relationship between TMB and prognosis and immune infiltration. It will increase the persuasiveness of the above results if we can further improve these aspects in the future.


In summary, lower TMB is correlated with worse survival outcomes and immune-related pathways, while higher TMB inhibits immune infiltration in SKCM patients. The four TMB-related immune gene model can effectively differentiate between high and low-risk patients, moreover, mutants of the four hub genes confer lower immune cell infiltration, which should be further validated.

Availability of data and materials

The TCGA-SKCM datasets:;

ImmPort Private Database:;

TIMER Database:



Area under the curve


Confidence interval


Differentially expressed genes


Disease Specific Survival


Gene ontology


Gene Set Enrichment Analysis


Hazard ratio


Immunocheck point blockade


Overall Survival


Progression Free Interval

ROC curve:

Receiver operating characteristic curve


Skin Cutaneous Melanoma


Single nucleotide polymorphism


Single nucleotide variant


The Cancer Genome Atlas


Tumor mutation burden


  1. 1.

    De La Cruz MM, Abdul Z, Shariff Z. The impact of a skin cancer diagnosis on travel insurance: a sun worshipper's dilemma. Clin Exp Dermatol. 2020.

  2. 2.

    Siegel JA, Yudkin JS, Craker K, Hwang A, Libby T. Uncapping the bottle: a proposal to allow full-sized sunscreens in carry-on luggage to promote sun protection and prevent skin cancer. J Am Acad Dermatol. 2020.

  3. 3.

    Khan AQ, Travers JB, Kemp MG. Roles of UVA radiation and DNA damage responses in melanoma pathogenesis. Environ Mol Mutagen. 2018;59(5):438–60.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  4. 4.

    Gupta R, Janostiak R, Wajapeyee N. Transcriptional regulators and alterations that drive melanoma initiation and progression. ONCOGENE. 2020;39(48):7093–105.

  5. 5.

    Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2018;68(6):394–424.

    Article  Google Scholar 

  6. 6.

    Siegel RL, Miller KD, Jemal A. Cancer statistics, 2019. CA Cancer J Clin. 2019;69(1):7–34.

    Article  Google Scholar 

  7. 7.

    Chen W, Zheng R, Baade PD, Zhang S, Zeng H, Bray F, et al. Cancer statistics in China, 2015. CA Cancer J Clin. 2016;66(2):115–32.

    Article  Google Scholar 

  8. 8.

    Luke JJ, Flaherty KT, Ribas A, Long GV. Targeted agents and immunotherapies: optimizing outcomes in melanoma. Nat Rev Clin Oncol. 2017;14(8):463–82.

    CAS  Article  PubMed  Google Scholar 

  9. 9.

    Schreuer M, Jansen Y, Planken S, Chevolet I, Seremet T, Kruse V, et al. Combination of dabrafenib plus trametinib for BRAF and MEK inhibitor pretreated patients with advanced BRAFV600-mutant melanoma: an open-label, single arm, dual-Centre, phase 2 clinical trial. Lancet Oncol. 2017;18(4):464–72.

    CAS  Article  PubMed  Google Scholar 

  10. 10.

    Long GV, Hauschild A, Santinami M, Atkinson V, Mandalà M, Chiarion-Sileni V, et al. Adjuvant Dabrafenib plus Trametinib in stage III BRAF-mutated melanoma. New Engl J Med. 2017;377(19):1813–23.

    CAS  Article  PubMed  Google Scholar 

  11. 11.

    Sanlorenzo M, Vujic I, Floris A, Novelli M, Gammaitoni L, Giraudo L, et al. BRAF and MEK inhibitors increase PD-1-positive melanoma cells leading to a potential lymphocyte-independent synergism with anti–PD-1 antibody. Clin Cancer Res. 2018;24(14):3377–85.

    CAS  Article  PubMed  Google Scholar 

  12. 12.

    Kunz M, Hölzel M. The impact of melanoma genetics on treatment response and resistance in clinical and experimental studies. Cancer Metast Rev. 2017;36(1):53–75.

    CAS  Article  Google Scholar 

  13. 13.

    Hu-Lieskovan S, Mok S, Homet Moreno B, Tsoi J, Robert L, Goedert L, et al. Improved antitumor activity of immunotherapy with BRAF and MEK inhibitors inBRAFV600E melanoma. Sci Transl Med. 2015;7(279):241r–79r.

    Article  Google Scholar 

  14. 14.

    Valpione S, Campana LG. Immunotherapy for advanced melanoma: future directions. Immunotherapy-UK. 2016;8(2):199–209.

    CAS  Article  Google Scholar 

  15. 15.

    Axelrod ML, Johnson DB, Balko JM. Emerging biomarkers for cancer immunotherapy in melanoma. Semin Cancer Biol. 2018;52:207–15.

    CAS  Article  PubMed  Google Scholar 

  16. 16.

    Dummer R, Ascierto PA, Nathan P, Robert C, Schadendorf D. Rationale for immune checkpoint inhibitors plus targeted therapy in metastatic melanoma: a review. Jama Oncol. 2020;6(12):1957.

    Article  Google Scholar 

  17. 17.

    Effern M, Glodde N, Braun M, Liebing J, Boll HN, Yong M, et al. Adoptive T cell therapy targeting different gene products reveals diverse and context-dependent immune evasion in melanoma. Immunity. 2020;53(3):564–80.

    CAS  Article  PubMed  Google Scholar 

  18. 18.

    Sha D, Jin Z, Budczies J, Kluck K, Stenzinger A, Sinicrope FA. Tumor Mutational Burden as a Predictive Biomarker in Solid Tumors. Cancer Discov. 2020;10(12):1808–25.

  19. 19.

    Jardim DL, Goodman A, de Melo Gagliato D, Kurzrock R. The challenges of tumor mutational burden as an immunotherapy biomarker. Cancer cell. 2021;39(2):154–73.

  20. 20.

    Heydt C, Rehker J, Pappesch R, Buhl T, Ball M, Siebolts U, Haak A, Lohneis P, Büttner R, Hillmer AM et al. Analysis of tumor mutational burden: correlation of five large gene panels with whole exome sequencing. Sci Rep-UK. 2020; 10(1). doi:

  21. 21.

    Ajona D, Ortiz-Espinosa S, Moreno H, Lozano T, Pajares MJ, Agorreta J, et al. A combined PD-1/C5a blockade synergistically protects against lung Cancer growth and metastasis. Cancer Discov. 2017;7(7):694–703.

    CAS  Article  PubMed  Google Scholar 

  22. 22.

    Razzak M. Anti-PD-1 approaches—important steps forward in metastatic melanoma. Nat Rev Clin Oncol. 2013;10(7):365.

    Article  PubMed  Google Scholar 

  23. 23.

    A Set of Transcriptomic Changes Is Associated with Anti–PD-1 Resistance. Cancer Discov. 2016; 6(5):471–472.doi:

  24. 24.

    Cristescu R, Mogg R, Ayers M, Albright A, Murphy E, Yearley J, et al. Pan-tumor genomic biomarkers for PD-1 checkpoint blockade–based immunotherapy. Science. 2018;362(6411):r3593.

    Article  Google Scholar 

  25. 25.

    Hugo W, Zaretsky JM, Sun L, Song C, Moreno BH, Hu-Lieskovan S, et al. Genomic and Transcriptomic features of response to anti-PD-1 therapy in metastatic melanoma. Cell. 2016;165(1):35–44.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  26. 26.

    Mayakonda A, Lin D, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 2018;28(11):1747–56.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  27. 27.

    Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47.

    Article  Google Scholar 

  28. 28.

    Yu G, Wang L, Han Y, He Q. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284–7.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  29. 29.

    Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci. 2005;102(43):15545–50.

    CAS  Article  PubMed  Google Scholar 

  30. 30.

    Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, Tamayo P. The molecular signatures database Hallmark gene set collection. Cell Syst. 2015;1(6):417–25.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  31. 31.

    Bhattacharya S, Dunn P, Thomas CG, Smith B, Schaefer H, Chen J, et al. ImmPort, toward repurposing of open access immunological assay data for translational and clinical research. SCI DATA. 2018;5(1):180015.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  32. 32.

    Li T, Fan J, Wang B, Traugh N, Chen Q, Liu JS, et al. TIMER: a web server for comprehensive analysis of tumor-infiltrating immune cells. Cancer Res. 2017;77(21):e108–10.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  33. 33.

    Newman AM, Steen CB, Liu CL, Gentles AJ, Chaudhuri AA, Scherer F, et al. Determining cell type abundance and expression from bulk tissues with digital cytometry. Nat Biotechnol. 2019;37(7):773–82.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  34. 34.

    Couzin-Frankel J. Cancer immunotherapy. Science. 2013;342(6165):1432–3.

    CAS  Article  PubMed  Google Scholar 

  35. 35.

    Van Allen EM, Miao D, Schilling B, Shukla SA, Blank C, Zimmer L, et al. Genomic correlates of response to CTLA-4 blockade in metastatic melanoma. Science. 2015;350(6257):207–11.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  36. 36.

    Howitt BE, Shukla SA, Sholl LM, Ritterhouse LL, Watkins JC, Rodig S, Stover E, Strickland KC, D Andrea AD, Wu CJ et al. Association of Polymerase e–mutated and microsatellite-instable endometrial cancers with Neoantigen load, number of tumor-infiltrating lymphocytes, and expression of PD-1 and PD-L1. Jama Oncol 2015; 1(9):1319, DOI:

  37. 37.

    Chan TA, Wolchok JD, Snyder A. Genetic basis for clinical response to CTLA-4 blockade in melanoma. New Engl J Med. 2015;373(20):1984.

    CAS  Article  PubMed  Google Scholar 

  38. 38.

    Rizvi NA, Hellmann MD, Snyder A, Kvistborg P, Makarov V, Havel JJ, et al. Mutational landscape determines sensitivity to PD-1 blockade in non–small cell lung cancer. Science. 2015;348(6230):124–8.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  39. 39.

    Hellmann MD, Ciuleanu T, Pluzanski A, Lee JS, Otterson GA, Audigier-Valette C, et al. Nivolumab plus Ipilimumab in lung Cancer with a high tumor mutational burden. New Engl J Med. 2018;378(22):2093–104.

    CAS  Article  PubMed  Google Scholar 

  40. 40.

    Rosenthal R, Cadieux EL, Salgado R, Bakir MA, Moore DA, Hiley CT, et al. Neoantigen-directed immune escape in lung cancer evolution. Nature. 2019;567(7749):479–85.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  41. 41.

    Desrichard A, Kuo F, Chowell D, Lee K, Riaz N, Wong RJ, et al. Tobacco smoking-associated alterations in the immune microenvironment of squamous cell carcinomas. J Natl Cancer Institute. 2018;110(12):1386–92.

    CAS  Article  Google Scholar 

  42. 42.

    Fredriksson NJ, Elliott K, Filges S, Van den Eynden J, Ståhlberg A, Larsson E. Recurrent promoter mutations in melanoma are defined by an extended context-specific mutational signature. PLoS Genet. 2017;13(5):e1006773.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  43. 43.

    Schreiber RD, Old LJ, Smyth MJ. Cancer Immunoediting: integrating Immunity's roles in Cancer suppression and promotion. SCIENCE. 2011;331(6024):1565–70.

    CAS  Article  PubMed  Google Scholar 

  44. 44.

    Rieckmann JC, Geiger R, Hornburg D, Wolf T, Kveler K, Jarrossay D, et al. Social network architecture of human immune cells unveiled by quantitative proteomics. Nat Immunol. 2017;18(5):583–93.

    CAS  Article  PubMed  Google Scholar 

  45. 45.

    Cao Y, Wang X, Jin T, Tian Y, Dai C, Widarma C, et al. Immune checkpoint molecules in natural killer cells as potential targets for cancer immunotherapy. Signal Transduction Targeted Therapy. 2020;5(1):250.

    Article  PubMed  Google Scholar 

  46. 46.

    Coleman E, Ko C, Dai F, Tomayko MM, Kluger H, Leventhal JS. Inflammatory eruptions associated with immune checkpoint inhibitor therapy: a single-institution retrospective analysis with stratification of reactions by toxicity and implications for management. J Am Acad Dermatol. 2019;80(4):990–7.

    CAS  Article  PubMed  Google Scholar 

  47. 47.

    Lynes J, Jackson S, Sanchez V, Dominah G, Wang X, Kuek A, et al. Cytokine microdialysis for real-time immune monitoring in glioblastoma patients undergoing checkpoint blockade. Neurosurgery. 2019;84(4):945–53.

    Article  PubMed  Google Scholar 

  48. 48.

    Routy B, Le Chatelier E, Derosa L, Duong CPM, Alou MT, Daillère R, et al. Gut microbiome influences efficacy of PD-1–based immunotherapy against epithelial tumors. Science. 2018;359(6371):91–7.

    CAS  Article  PubMed  Google Scholar 

  49. 49.

    Schumacher TN, Schreiber RD. Neoantigens in cancer immunotherapy. Science. 2015;348(6230):69–74.

    CAS  Article  Google Scholar 

Download references


Not applicable.


This study has supported by Health Commission Scientific Research Program of Chongqing (NO:2012–2-056). The funding agency had no role in the study design, collection, analysis, or interpretation of data, in the writing of this manuscript, or in the decision to submit for publication.

Author information




(I) Conception and design: LWW and ZJY; (II) Performed data analysis and created the figures/tables: FC and RL; (III) Data analysis and interpretation: LWW, LS, and GSZ; (IV) Manuscript drafting and editing: All authors; (V) Final approval of manuscript: All authors.

Corresponding author

Correspondence to Zhengjian Yan.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Figure S1.

Workflow of this study.

Additional file 2: Figure S2.

The genecloud plot and the coincidences&exclusive relationship plot among the mutated genes.

Additional file 3: Figure S3.

Venn analysis of 1812 immune genes and 504 DEGs.

Additional file 4: Figure S4.

Correlation and distribution analysis of the infiltration degree of 22 immune cells in the high-TMB and low-TMB groups.

Additional file 5: Figure S5.

K-M analysis based on the infiltration levels.

Additional file 6: Table S1-S7.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Wang, L., Chen, F., Liu, R. et al. Gene expression and immune infiltration in melanoma patients with different mutation burden. BMC Cancer 21, 379 (2021).

Download citation


  • TCGA
  • SKCM
  • TMB
  • ICB
  • Survival analysis