Ronghua Yang1†
Zhengguang Wang
Sitong Zhou- 1Department of Burn Surgery and Skin Regeneration, The First People’s Hospital of Foshan, Foshan, China
- 2Department of Orthopedics, The First Affiliated Hospital of China Medical University, Shenyang, China
- 3Department of Dermatology, The First People’s Hospital of Foshan, Foshan, China
- 4Department of Anesthesiology, The First People’s Hospital of Foshan, Foshan, China
- 5Department of Burns, Nanfang Hospital, Southern Medical University, Guangzhou, China
- 6Department of Oncology, Affiliated Hospital of Chengdu University, Chengdu, China
Skin cutaneous melanoma (SKCM) is a highly aggressive and resistant cancer with immense metabolic heterogeneity. Here, we performed a comprehensive examination of the diverse metabolic signatures of SKCM based on non-negative matrix factorization (NMF) categorization, clustering SKCM into three distinct metabolic subtypes (C1, C2, and C3). Next, we evaluated the metadata sets of the metabolic signatures, prognostic values, transcriptomic features, tumor microenvironment signatures, immune infiltration, clinical features, drug sensitivity, and immunotherapy response of the subtypes and compared them with those of prior publications for classification. Subtype C1 was associated with high metabolic activity, low immune scores, and poor prognosis. Subtype C2 displayed low metabolic activity, high immune infiltration, high stromal score, and high expression of immune checkpoints, demonstrating the drug sensitivity to PD-1 inhibitors. The C3 subtype manifested moderate metabolic activity, high enrichment in carcinogenesis-relevant pathways, high levels of CpG island methylator phenotype (CIMP), and poor prognosis. Eventually, a 90-gene classifier was produced to implement the SKCM taxonomy and execute a consistency test in different cohorts to validate its reliability. Preliminary validation was performed to ascertain the role of SLC7A4 in SKCM. These results indicated that the 90-gene signature can be replicated to stably identify the metabolic classification of SKCM. In this study, a novel SKCM classification approach based on metabolic gene expression profiles was established to further understand the metabolic diversity of SKCM and provide guidance on precisely targeted therapy to patients with the disease.
Introduction
Skin cutaneous melanoma (SKCM) is the deadliest type of skin cancer due to its high metabolic and metastatic rates, accounting for more than 80% of skin cancer-related deaths (Guy et al., 2015; Bolick and Geller, 2021). In the past few years, the advancement of immune checkpoint blockade agents has become a pillar of SKCM therapy, remarkably boosting the therapeutic outcomes. However, the response of patients with SKCM to immunotherapy is heterogeneous, with approximately 50% of them experiencing unfavorable responses (Robert et al., 2011; Schadendorf et al., 2015). Therefore, it is imperative to uncover the latent molecular mechanisms of SKCM heterogeneity to develop precise immunotherapies and determine the populations that would benefit the most from them.
Skin cutaneous melanoma is characterized by prominent metabolic plasticity, a feature that results from the activation of oncogenic pathways due to a high frequency of somatic mutations (Alexandrov et al., 2020). SKCM has been classified into four genomic subtypes, BRAF subtype, RAS subtype, NF1 subtype, and triple wild type based on the somatic mutations in these genes and their ratios (2015). These intrinsic oncogenes contribute to the metabolic conversion of SKCM, which leads to a high degree of plasticity and adaptation of melanoma to unfavorable conditions (Ratnikov et al., 2017). Moreover, the transformed metabolic microenvironment can reprogram the function of immune cell subpopulations, allowing melanoma to evade the immune system (Bristot et al., 2020).
A recent proteogenomic research divided patients with SKCM into two subgroups according to whether they responded to immunotherapy against PD-1 or TIL (Harel et al., 2019), namely, the responder and the non-responder subgroups. Differential protein expression analysis found that differences in mitochondrial metabolism were responsible for the response of patients to immunotherapy rather than the treatment protocol. In summary, the heterogeneity of melanoma metabolism is an important reason for the poor efficacy of immunotherapy. Therefore, this study was conducted to categorize SKCM from a metabolic viewpoint to reveal its heterogeneity.
In the present study, we performed a systematic examination of the diverse metabolic signatures of SKCM using a screened metabolic gene based on a non-negative matrix factorization (NMF) clustering algorithm and identified three distinct metabolic subtypes. In this process, The Cancer Genome Atlas (TCGA)-SKCM cohorts were merged into a metadata set for clustering based on the expression of metabolism genes. Additional processed microarray profiles of Gene Expression Omnibus (GEO) SKCM samples were used for external validation. Unsupervised transcriptomic analysis identified three subtypes of SKCM, namely, C1, C2, and C3. In addition, by comparing transcriptomic data from patients with different metabolic subtypes, differentially expressed genes (DEGs) were retrieved. We estimated the prognostic difference, transcriptome features, relationships with metabolic signatures, tumor microenvironment features, immune infiltration, clinical traits, somatic mutation signatures, immunotherapy, and drug sensitivity of the SKCM subtypes, and a comparison was made with previously established classifications. Finally, a 90-gene classifier was used to determine SKCM classification. This research may also provide in-depth insights into tumor–immune cell interactions, showing considerable promise for the clinical therapeutic interventions of patients with SKCM.
Materials and Methods
Patients and Samples
Gene expression profiles of SKCM, including TCGA-SKCM (Cancer Genome Atlas Network, 2015), GSE54467, and GSE65904, were obtained from three independent cohorts of patients. In addition, only SKCM samples were retained for further analysis. TCGA-SKCM project was downloaded using the TCGAbiolinks package (Colaprico et al., 2016) and converted to the TPM format for subsequent analysis. The annotation of Ensembl ID for protein-coding mRNAs was transformed to the gene symbol based on the GENCODE gene model (GENCOED27). Then the batch effect from the different datasets was corrected using the ComBat package (Zhang et al., 2020) in the SKCM cohorts. Clinical data regarding the disease, including age, sex, tumor stage, and survival information, were retrieved from TCGA Pan-Cancer Clinical Data Resource (TCGA-CDR), and the clinical characteristics of the TCGA-SKCM patients are shown in Supplementary Table 1, of which only the overall survival (OS) information was obtained for further data processing. In addition, the copy number mutation data of TCGA-SKCM cohorts were downloaded from the GDAC-Firebrowse website1.
Identification of Skin Cutaneous Melanoma Molecular Subtypes by Non-negative Matrix Factorization Clustering
The NMF clustering (Possemato et al., 2011) algorithm was used to cluster the SKCM samples. The 2,752 metabolism-related genes that encode all the well-known human metabolic enzymes and transporters were selected for follow-up screening. First, the metabolism-related genes that were significantly correlated with OS time were subjected to Cox survival regression using the survival R package. Then unsupervised NMF clustering (Gaujoux and Seoighe, 2010) was performed based on the TCGA-SKCM cohort, and validation was performed from the integrated cohorts from GSE54467 and GSE65904 using the same selected candidate genes. K values were chosen where the magnitude of the cophenetic correlation coefficient started to decrease with the optimum number of clusters (Brunet et al., 2004). Next, we evaluated the similarity of subtype classification between independent cohorts based on the expression profiles of mRNAs by employing the class mapping analysis (SubMap) (Gene pattern) method to assess whether the subtypes were analyzed in the training, and validation sets were significantly correlated. Simultaneously, the mRNA expression data of the abovementioned candidate genes were analyzed to verify the subtype distributions using the T-distributed stochastic neighbor embedding (t-SNE) method.
Gene Set Variation Analysis
Gene set variation analysis (GSVA) (Hänzelmann et al., 2013) is a gene set enrichment method that computes an estimated fraction of certain pathways or signatures of different clusters based on expression profiles. The relevant metabolism (Rosario et al., 2018) and carcinogenesis- (Sanchez-Vega et al., 2018)-relevant pathway gene sets were obtained from previous studies, and the GSVA R package was used to investigate the gene set differences between samples. Subsequently, metabolism gene scores were obtained for differential analysis using the Limma package (Ritchie et al., 2015) in R software, and differentially expressed signatures were screened out with the following threshold (| log2FC| > 0.2, adjusted p < 0.05).
Estimation of the Tumor Microenvironment Signatures
The microenvironment cell population counter (MCPcounter) (Becht et al., 2016) was used to estimate the number of infiltrated immune cell populations and two non-immune stromal cell populations (immune cell types: T cells, CD8 + T cells, natural killer cells, cytotoxic lymphocytes, B-cell lineage, monocytic lineage, myeloid dendritic cells, and neutrophils; stromal cell types: endothelial cells and fibroblasts). Furthermore, another approach applied to quantify tumor immune components was the single-sample gene set enrichment analysis (ssGSEA) method (Barbie et al., 2009), which calculates enrichment scores representing the degree to which genes in a particular gene set are coordinately upregulated or downregulated within a single sample. In particular, six immune cell populations, including regulatory T cells (Treg), helper T cells 1 (Th1), helper T cells 2 (Th2), helper T cells 17 (Th17), central memory T cells, and effective memory T cells (Tem), were analyzed using the GSVA R package. In addition, the ESTIMATE algorithm (Yoshihara et al., 2013) was applied to computer immune and stromal scores in different subtypes, thereby, reflecting the features of the tumor microenvironment.
Characterization of Skin Cutaneous Melanoma Subtypes
After data normalization, differentially expressed genes (DEGs) between different SKCM subtypes were identified using the Limma package (| log2FC| > 1 and p < 0.01). The gene signature set files “c2.cp.kegg.v6.2. symbols.gmt” and “h,all,v60.2.symbols” were downloaded from the Molecular Signature Database (MSigDB)2. GSEA was then applied to investigate the pathway and functional enrichment using the Clusterprofiler R package (Yu et al., 2012) with the significance threshold set to an adjusted p < 0.05. Furthermore, previously published molecular classifications of SKCM were predicted using the nearest template prediction (NTP) analysis from gene pattern modules, and then the prediction outcome was compared with our classification.
Construction and Validation of the Skin Cutaneous Melanoma Gene Classifier
To identify specific genes in the SKCM subtypes, we screened genes with statistically significant differences in the different subclasses according to the following criteria: for adjusted p < 0.01 and absolute log2 FC > 2. Only the genes with a significantly different expression in all three possible parameters were considered as subclass-specific genes. The top 30 genes with the largest log2FC values in each subtype (only genes with log2FC > 0 were selected) were further used to generate the prediction models such that a 90-gene subtype classifier was created. Next, we used the NTP algorithm to predict the subclasses of the 90-gene signature in GSE14520 and compared them with the previous classification results derived from the NMF algorithm.
Prediction of the Efficacy of Each Subtype of Immunotherapy and Targeted Therapy
We used data from patients with melanoma treated with immunotherapy to indirectly predict the efficacy of immunotherapy in melanoma subclasses by measuring the similarity of gene expression profiles between the subclasses determined in this study and those in patients with melanoma based on SubMap analysis (gene pattern). Furthermore, we downloaded and performed SubMap analysis from the genomics of drug sensitivity in cancer (GDSC) database (Yang et al., 2013) to investigate its drug sensitivity.
Gene Ontology and KEGG Analyses
We performed Gene Ontology (GO) and KEGG enrichment analysis for the different subclasses of differentially expressed genes, where the GO includes biological process (BP), molecular function (MF), and cellular component (CC). GO analysis of differentially expressed genes was conducted using DAVID (Huang et al., 2007; Huang da et al., 2009) (FDR < 0.1) and visualized using the ggplot23 R package and Goplot package (Walter et al., 2015).
We performed the KEGG pathway enrichment analysis using the KOBAS 3.0 database4 (Wu et al., 2006) for the integration of DEGs with Ensemble ID in the differential gene list and then obtained pathway enrichment lists where pathways with p < 0.05 were considered significantly enriched.
Mutation Analysis Using a 90-Gene Classifier
The gene mutation and gene copy number data of the 90-gene classifier in 32 TCGA pan-cancer databases were retrieved from the cBioportal (Gao et al., 2013) web portal. The mutation status of these genes in these databases were analyzed, and a bar chart showing the distribution ratio of the 90 genes in TCGA pan-cancer databases by mutation type, fusion, amplification, deep deletion, and multiple alterations were constructed.
Transcription Factor Prediction
Transcription factors often modulate several metabolic genes that are closely functionally associated. Therefore, we utilized the NetworkAnalyst network platform (Zhou et al., 2019) to analyze and predict the transcription factors that are most likely to regulate the 90 genes included in the classifier and construct a molecular interaction network. The targeted gene–transcription factor interaction network also contains the TF–mRNA–miRNA molecular regulatory network data obtained from the RegNetwork information library.
Expression of Solute Carrier Family 7 Member 4 in Skin Cutaneous Melanoma and Kaplan–Meier Analysis
RNA-seq data were downloaded from UCSC XENA in TPM format, and GTEx was downloaded from UCSC XENA5 and processed by the Toil project (Vivian et al., 2017). The Wilcoxon rank-sum test was used to compare solute carrier family 7 member 4 (SLC7A4) expression in normal TCGA and GTEx skin tissues and tumor samples from SKCM in TCGA. We categorized patients into high and low SLC7A4 expression groups according to the median value of their SLC7A4 expression (TPM format), following which Kaplan–Meier analysis of OS was performed for each group and visualized using the R package “survminer.”
Patient Tissue Specimen Collection and Immunohistochemistry Validation
Melanoma tissues were collected from 25 melanoma patients of Han Chinese ethnicity from 2015 to 2020. Informed consent was obtained from each patient, and the study was approved by The Foshan Subject Review Board of the First People’s Hospital. Paraffin-embedded tissues were sectioned at 4-mm thickness for immunohistochemistry (IHC) analysis. Antigen retrieval was performed by incubating the samples in citrate buffer (pH 6.0) for 15 min. After blocking with a mixture of methanol and 0.75% hydrogen peroxide, sections were incubated overnight with a primary antibody (SLC7A4, Proteintech, 1:200) followed by incubation with a secondary antibody conjugated with horseradish peroxidase (HRP; goat anti-rabbit IgG, 1:500, Cell Signaling Technology). The sections were washed three times with PBS and incubated with AEC (ZSGB-BIO). Staining was performed as described previously (Zhou et al., 2021). Tissues were examined through the cross-product (H score) of the percentage of tumor cell staining at each of the three staining intensities, and the staining score was graded through the H score as follows: low, H score = 0–100; moderate, H score = 101–200; and high, H score = 201–300.
Statistical Analysis
All data processing and analyses were performed using Excel (Microsoft) and R software (version 4.0.2). Unpaired Student’s t-test was employed for the comparison of two groups with non-normal distribution, and the Mann–Whitney U-test was used for the comparison of two groups with non-normal distribution. One-way analysis of variance (ANOVA) and Kruskal–Wallis tests were performed for comparisons among three groups. Contingency table (χ2) variables used the χ2 test for statistical significance. Survival analysis was performed using the Kaplan–Meier method, and the log-rank test was used for comparison. Univariate Cox proportional risk regression models were used to evaluate risk ratios for univariate analyses. Two-tailed p < 0.05 was considered statistically significant.
Results
Non-negative Matrix Factorization Determined Three Subtypes of Skin Cutaneous Melanoma
Before conducting the NMF algorithm analysis of SKCM, we first utilized the ComBat algorithm to remove the batch effect for different SKCM cohorts and chart the PCA after batch effect removal (Figure 1A). A total of 2,752 previously reported metabolism-related genes were selected as the basis for NMF analysis. We then adopted univariate Cox regression for metabolism-related genes in the metadata set to identify prognostic genes associated with OS (p < 0.1), and a total of 517 candidate genes were obtained. We then extracted TCGA data expression profiles of SKCM, clustering the 517 candidate genes using the NMF clustering algorithm, and the non-negative matrix decomposition (NMF) with two to six clusters was plotted (Figure 1B). The best clustering number k value was established by computing the clustering correlation coefficient, and k = 3 was regarded as the optimal clustering number. The consensus NMF was performed again with several decompositions of 3, defining three subtypes C1 (n = 113), C2 (n = 103), and C3 (n = 234). To validate the clustering subtypes, we performed t-SNE dimension reduction to reduce the feature dimensionality for all samples of metabolic genes in the expression profile of the test TCGA dataset and showed the different two-dimensional distribution pattern plots for the three types of samples (Figure 1C), and the clinical characteristics of TCGA-SKCM are shown in Supplementary Table 1. We also found that our subtypes were largely consistent with the two-dimensional t-SNE distribution pattern. The same consensus NMF for the same set of metabolic genes was performed for the consolidated validation set (GSE54467 and GSE65904), and three classes C1 (n = 93), C2 (n = 98), and C3 (n = 98) were obtained. Finally, the subclass algorithm was used to identify the subtype matching model of TCGA and validation sets, and it was decided that TCGA-C1 = GEO-C1, TCGA-C2 = GEO-C3, and TCGA-C3 = GEO-C2. We also utilized survival information from the three cohorts to conduct a subtype survival analysis of the SKCM subset. In TCGA-SKCM cohort, the results showed a significantly higher OS in C2 than in C1 and C3 (log-rank test p < 0.0001, Figure 1D), and the same survival differences were also verified in the GEO validation dataset (GSE65904 and GSE54467). Similar survival outcomes were observed; however, the prognosis of C2 in the validation set worsened, and the median survival time was more similar to that of C1 (Figure 1E).
Figure 1. Identification of skin cutaneous melanoma (SKCM) subtypes using non-negative matrix factorization (NMF) consensus clustering in the metadata set. (A) Visualization of principal component analysis of combined expression profiles of three cohort data sets. (B) Phenotype correlation coefficients for NMF clustering of 816 metabolism-related genes at k = 2–5. (C) t-SNE analysis rendered in support of classification into three SKCM subtypes. (D) Overall survival analysis of three subtypes (C1, C2, and C3) in two independent cohorts—test set (TCGA-SKCM) and (E) Validation set (GSE65904, GSE54467). Statistical significance of differences was determined using the log-rank test.
Association of Skin Cutaneous Melanoma Subtypes With Metabolism-Related Signatures
Considering that the SKCM subclass classification was based on metabolism-related genes, we further investigated whether different metabolic signatures are present in the different subclasses. First, we scored the metabolic pathways (gene sets were acquired from the reported paper) using the GSVA R package in the subtype cohorts. In TCGA-SKCM cohort, GSVA enrichment scores of the metabolic pathways were estimated, and a cross-group Limma difference test was performed using two groups of “subtype n vs. other subtypes” to confirm subtype-specific differential metabolic pathways; the screening standard for GSVA enrichment was | logFC| > 0.2, adjusted p < 0.05, and heat maps were constructed for visualization (Figure 2A).
Figure 2. Correlation between metabolically relevant features and immune infiltration in SKCM subtypes. (A) Heat map of specific metabolism-related signatures in three subtypes. (B) The heat map depicted the richness of immune and stromal cell populations in C1, C2, and C3 by MCPcounter and ssGSEA.
The heat map showed that C1, C2, and C3 had 21, 3, and specific metabolic signatures, respectively. Among them, C1 demonstrated a distinct metabolic signature, with 21 metabolic pathways significantly upregulated. Both C2 and C3 were significantly downregulated in these pathways. These results demonstrated that each subtype was enriched in unique metabolic pathways and had dissimilar metabolic levels.
We also counted the GSVA enrichment points of carcinogenesis-related pathways and plotted a box line among groups. To further investigate their subtype characteristics, we selected a collection of 11 carcinogenesis-related pathways and quantified them using the GSVA algorithm. The results showed significant between-group differences among the three subgroups in cell cycle, HIPPO, MYC, NRF2, PI3K, TGF-β, TP53, WNT, and angiogenesis, which are carcinogenic signaling pathways (Figure 3A), indicating a close connection between our subtypes and carcinogenesis. The results revealed that C1 had a significantly strong cell cycle and WNT signature than C2 and C3, C2 displayed increased expression of components of the angiogenesis pathway, and C3 was particularly enriched in the HIPPO, NRF2, PI3K, TGF-β, and TP53 pathways.
Figure 3. Tumor microenvironment signatures of three subtypes in the TCGA-SKCM. (A) Box plots of characteristic marks discriminating different subtypes of SKCM progression-related features. Box-line plots of the immune score (B) and stromal score (C) estimate three subtypes. For the box-line plot, the lines in the box represent the median, the bottom, and top of the box, which are the 25th and 75th percentiles (interquartile range), and the vertical lines represent 1.5 times the interquartile range. Statistical differences were compared by the Kruskal–Wallis test, and p-values are marked with an asterisk above each box line plot (ns indicates no significance, **p < 0.01, ****p < 0.0001). (D) The boxplot showed the accumulation of immune and stromal cell populations discriminated by different subtypes. (E) Expression level (normalized transcripts per million) of 14 immune checkpoint genes in three SKCM subtypes. The difference was verified statistically through the Kruskal–Wallis test, and the p-values are noted with asterisks at the top of each boxplot (ns stands for no significance, *p < 0.05, ***p < 0.001).
Correlation Between Skin Cutaneous Melanoma Subtypes and Immune Infiltration
As significant differences in immune scores were detected between SKCM subtypes, we surveyed immune cell infiltration in SKCM subtypes to assess their immunological landscape. We utilized the MCP-counter and ssGSEA algorithms to compute the abundance of 16 immune infiltration cells and presented them in an immune heat map (Figure 2B). We also plotted box plots of inter-group differences in immune cells, and the results displayed significant inter-group differences in all types of immune cell populations among the three subtypes (Figure 3D). Notably, the box plots revealed that C2 showed a significantly higher enrichment immune score than C1 and C3 in almost all immune cells, except for Th17 cells, in which both C1 and C2 had significantly higher Th17 cell enrichment scores than C3. Of these results, C2 was enriched with more immune cells, which is consistent with the finding that C2 had the highest immune score among the three subtypes. With the current widespread use of immune checkpoint inhibitors (ICIs) in clinical trials and for the treatment of advanced melanoma (Woods et al., 2016), we explored the correlation between the expression of 13 classically targeted immune checkpoint genes in these subtypes, which are currently based on immunotherapy inhibitors in clinical trials or licensed for certain cancer types. The results showed significant differences between groups, while C2 displayed a higher expression of nine immune checkpoint genes than C1 and C3, except for FGL1, CTAG1B, and MAGEA4 (Figure 3E).
Relevance of Skin Cutaneous Melanoma Subtypes to the Clinical Characteristics of Patients and in TCGA and GEO Datasets
To investigate the associations between subtypes and clinical characteristics, we analyzed the clinical tumor pathology variables associated with subtypes based on TCGA-SKCM (Figure 4A) and GEO validation set (Figure 4B) cohorts to construct the clinical information heat map of subtypes. The results showed an independence test for discrepancies between clinicopathological characteristics and metabolic subtypes. Furthermore, we matched this metabolic classification with the previously reported subtypes of SKCM in the literature, including mutation subtypes (BRAF, NF1, RAS, and Triple-WT), MethTypes (CpG island-methylated, Hyper-methylated, Hypo-methylated, and Normal-like), the MIR cluster (four types), the Protein Cluster (four classes), and the classification of TCGA (Immune, Keratin, and MITF-low) (Figure 4A and Supplementary Tables 2, 3). In TCGA-SKCM, the C1 subtype was significantly correlated with protein cluster 2 (p = 0.006), MIR type 3 (p < 0.001), hypomethylated type (p < 0.001), and keratin type (p < 0.001). The C2 subtypes associated with protein cluster 1 (p = 0.006), MIR type 2 (p < 0.001), normal-like methylated type (p < 0.001), and immune type (p < 0.001). The C3 subtype was related to protein cluster 3 (p = 0.006), MIR type 3 (p < 0.001), hyper-methylated type (p < 0.001), and immune type (p < 0.001).
Figure 4. Clinical characteristics of three subtypes in the TCGA-SKCM and GSE14520 cohorts. (A) Correlation of C1, C2, and C3 subtypes with clinical features and previous SKCM subtypes in the TCGA cohort. (B) Correlation of our classification with clinical features in the validation GSE14520 set cohorts.
Association of Skin Cutaneous Melanoma Subtypes With the Tumor Microenvironment
To further evaluate whether the subtypes were associated with the tumor microenvironment, we estimated the immune and stromal scores using the ESTIMATE algorithm for each group and constructed box-violin plots (Figures 3B,C). The results showed statistically significant differences in immune scores between the three groups, while the stromal scores were not statistically significant among the three groups. C2 had the highest immune score among the three groups.
The melanoma mutation landscape has been demonstrated to lead to alterations in the tumor microenvironment (Nassar and Tan, 2020) and immunotherapeutic response (Havel et al., 2019). Next, we investigated whether somatic mutation frequencies varied across SKCM subtypes and examined the mutation patterns in these subtypes. We applied the maftools package (Mayakonda et al., 2018) to estimate TCGA-SKCM-driven gene mutations and mapped waterfall plots of clusters within each group. High frequencies of mutations were observed for BRAF, COL5A1, NRAS, MECOM, NF1, ARID2, TP53, and CDKN2A in SKCM subtypes (Figure 5A and Supplementary Table 4). We then classified the total number of mutations and the predicted neoantigens (Figure 5B). The tumor mutation burden (TMB) was calculated for each metabolic subtype (Figure 5C), and box-line violin plots were constructed separately. Subsequently, we determined the frequency of amplification (Figure 5D) and deletion (Figure 5E) for the three subtypes using online GISTIC 2.0 analysis (Mermel et al., 2011) to construct box-line violin plots of amplifications vs. deletions for each subtype.
Figure 5. Relationship of SKCM subtypes with tumor mutation burden characteristics. (A) Driver-type oncogenic mutations based on TCGA-SKCM typing with intra-group aggregation waterfall plots. (B–E) Violin plots of the number of mutation-predicted neoantigens (B,C) and copy number aberrations (D,E) box lines for SKCM subtypes. Wilcoxon rank-sum test was used to compare statistical differences (ns indicates not significant).
Finally, we further reviewed the chromosomal segment values of the three subtypes to ascertain whether there were significant copy number alterations by performing online GISTIC2.0 analysis on genePattern6 to map the copy number change cytoband of each subtype. Cytoband revealed an overall description of the copy number variation in each subgroup (red representing gains and blue representing losses) (Figures 6A–C). As expected, our profiling demonstrated that significant copy number alterations emerged in the three subtypes, including those observed in the chromosomal region of 9p21.3 (CDKN2A) (Ghiorzo et al., 2006); amplification at 11q13.3 (CCND1) (Gibcus et al., 2007) in C1 and C2; amplification at 22q13.2 (TOB2) (Thanasai et al., 2006) in G1 and G3, in which TOB2 was significantly amplified in C1; and a major amplification at 5p13.3 (TERT) (Peifer et al., 2015) in the C1 subtype. While the cytobands for C1 and C3 displayed more regional amplifications and deletions than C2, this could also be an essential explanation for the superior prognosis of C2 over the C1 and C3 subtypes. Thus, the alteration of copy number may be a dominant mechanism responsible for the differences in metabolism and prognosis among the three subtypes.
Figure 6. The landscape of somatic copy-number alterations in the three subtypes. (A–C) Cytoband showed the genomic copy-number change from three subtypes. Amplification (red) and deletion (blue) of each regional peak are shown and highlighted.
Ninety-Gene Classifier and Performance Validation
To build a subtype classifier for clinical application, we further selected subtype-related signature genes. The results of the Limma differential expression analysis were based on the whole expression map “subtype n vs. other subtypes.” Using MOVICS (Lu et al., 2020) package analysis, we retrieved the top 30 genes specifically upregulated for each metabolic subtype as biomarkers and then constituted clinical models and plotted correlation heat maps (Figure 7A and Supplementary Table 4). Consequently, we derived a 90-gene classifier and executed a consistency test by operating the NTP algorithm to predict the metabolic subtype attribution of each sample in the TCGA cohort as well as the GEO test set cohort, and plotted the heat map of the true subtype if it matched the predicted subtype (Figures 7B,C). The results were largely consistent between NMF and NTP in the three different subtypes of the test and validation sets (κ = 0.631 p < 0.001, κ = 0.714 p < 0.001), indicating that this 90-gene signature could be replicated to identify the metabolic classification of SKCM.
Figure 7. Identification of predictive metabolism-gene classifier and prediction immunotherapeutic response. (A) The heat map exhibited the expression level of the 90-gene classifier in C1, C2, and C3. (B) Consistency of SKCM metabolism subtype predictions between the 90-gene classifier and the original NMF-based predictions in TCGA-SKCM. (C) Consistency of SKCM metabolism subtype predictions between the 90-gene classifier and the original NMF-based predictions in GSE54467 and GSE65904. (D) C2 is likely sensitive to the PD-1 receptor inhibitor (nominal p = 0.01) by SubMap analysis. (E) Top 12 box plots of predicted IC50 values based on GDSC database drugs in three subtypes of TCGA-SKCM dataset.
Specific Sensitivity of Skin Cutaneous Melanoma Subtypes to Immunotherapy and Potentially Targeted Therapies
On the one hand, the different patterns of immune infiltration and expression levels of immune checkpoint genes in different SKCM subtypes suggest the need to further investigate the possibility of an immunotherapeutic response. Therefore, we matched the expression profiles of the three subtypes using a subclass algorithm to ascertain the degree of similarity of TCGA metabolic subtypes to the response profiles of 47 patients with melanoma receiving immunotherapy (Hoshida et al., 2007; Roh et al., 2017). The results indicate that the TCGA-C2 subtype was more likely to be responsive to immunotherapy with anti-PD1 treatment (Figure 7E) (Bonferroni correction, p = 0.01).
On the other hand, to identify potential anti-melanoma drugs that are associated with ICIs, we attempted to find potentially sensitive and selective chemotherapy drugs using the GDSC drug sensitivity database. We compared drug sensitivity for more than 100 drugs in the GDSC database, and the top 12 differential response drugs were plotted and listed according to the Kruskal-Wallis test (Figure 7D). We detected a significant difference in the estimated IC50 values between the three subtypes and found that C2 may be more sensitive to chemotherapeutic drugs.
Functional and Pathway Enrichment Analysis of the 90-Gene Classifier
Gene Ontology and KEGG enrichment analyses were conducted for the three group DEGs using the DAVID database, and GO enrichment analyses were performed for each of the three GO categories: BP, MF, and CC. The ggplot2 package was used for visualization, and the top 10 enrichment results are shown (Supplementary Figure 1 and Supplementary Table 6–8). The GO enrichment results of the three groups of DEGs were mainly enriched in the metabolism process, immune-specific BP, protein binding, cell cycle, DNA damage and repair, and cancer-associated biological processes. By taking the intersection of the enrichment analysis results of the three groups of differential genes in the KOBAS database, we identified six KEGG pathways related to oncogenesis development and metabolism for presentation. The red dots represent the upregulated enriched genes, mainly enriched in apoptosis, cell cycle, HIF-1 signaling pathway, human T-cell leukemia virus 1, cancer infection pathway, PI3K–Akt signaling pathway, and carcinogenesis biological pathway (Supplementary Figure 2 and Supplementary Table 9).
Mutation Analysis of the Gene Classifier in TCGA Pan-Cancer Database
Based on the gene mutation and gene copy number data of the 90-gene classifier in 32 TCGA pan-cancer databases retrieved from the cBioportal web portal, the mutation rates of the genes in the classifier were analyzed, and the distribution ratio of each mutation type of the gene classifier in each pan-cancer tumor, including gene mutation, fusion, amplification, deep deletion, and multiple alterations, is shown in a bar chart (Supplementary Figure 3A). Among them, the mutation rate exceeded 50% in most tumor cancer types in TCGA, and the mutation rate was higher in lung squamous cell carcinoma, esophageal adenocarcinoma, stomach adenocarcinoma, SKCM, and bladder urothelial carcinoma, with an alteration frequency of more than 80%.
The TF–mRNA–miRNA Network Construction
Many closely related genes are often subject to simultaneous regulation by specific transcription factors; therefore, we utilized the NetworkAnalyst network platform to predict the transcription factors most likely to regulate the 90 genes in the classifier by analyzing them and constructing a TF–mRNA–miRNA molecular interaction regulatory network (Supplementary Figure 3B).
Expression Signature, Prognostic Value, and Preliminary Experimental Validation of Solute Carrier Family 7 Member 4
The 90-gene classifier distribution in the three subtypes and the fold changes are plotted in the heat map shown in Figure 7A and Supplementary Table 5. Of these, SLC7A4 was not reported to be associated with melanoma among the 90 genes. Therefore, we performed an analysis and preliminary experimental validation of the role of SLC7A4 in SKCM. SLC7A4 was overexpressed in SKCM than in the adjacent normal skin tissues (Figure 8A). The OS was poorer in patients with high SLC7A4 expression than in those with low SLC7A4 expression (hazard ratio [HR] = 1.34 [1.02–1.76], p = 0.038, Figure 8B). SLC7A4 expression was determined in tissue samples using IHC staining. The results showed that melanoma tissues presented the highest expression of SLC7A4, followed by benign nevi, whereas normal skin tissues showed the lowest expression (Figures 8C,D). Generally, the intensity of SLC7A4 staining increased with the progression of melanoma. These results demonstrated that SLC7A4 plays a hub role in SKCM progression and, to some extent, proved the accuracy of this gene classifier.
Figure 8. The expression signature, prognostic value, and preliminary validation of solute carrier family 7 member 4 (SLC7A4). (A) The expression of SLC7A4 was higher in SKCM than in normal skin tissues. (B) The overall survival analysis of SLC7A4 in SKCM. (C,D) IHC staining demonstrated that the expressions of SLC7A4 staining increased with the progression of the disease. Normal skin tissues (n = 20) showed the lowest expression, followed by benign nevi (n = 20), and malignant melanoma tissues (n = 25) showing the highest expression. IHC, immunohistochemistry. IHC stain, AEC, original magnification × 100, scale bar = 0.1 mm (inset, IHC stain, AEC, original magnification ×400).
Discussion
Metabolic reprogramming is a key hallmark of cancer (Pavlova and Thompson, 2016). On the one hand, SKCM is characterized by its marked metabolic plasticity, and its development is associated with an important relationship between carcinogenesis and energy metabolism (Ruocco et al., 2019). On the other hand, metabolic reprogramming of melanoma alters the subset function of immune cells, enabling melanoma to evade the immune system. Thus, there is an urgent need to better understand the metabolic signatures of tumor heterogeneity in patients with melanoma to provide precise immunotherapy selection strategies and characterize the mechanisms of drug resistance. Herein, we performed a thorough classification of the metabolic profiles of SKCM specimens. Our findings revealed that SKCM could be categorized into three different metabolically relevant subtypes, and we verified the reproducibility of these subtypes in the GEO test set. Each subtype was characterized by different metabolic signatures, prognosis values, clinical parameters, tumor microenvironment features, immune cell infiltration, function and pathway enrichment, somatic copy number alterations, response to immunotherapy, and drug sensitivity. The C1 subtype displayed high metabolism levels, resembling the keratin subtype, with an advanced pathological stage and low immune infiltration. Furthermore, C1 had the worst prognosis and the shortest DFS among SKCM patients. In contrast, C2 was abundant in immune signatures with a relatively high expression of immune checkpoint genes and high immune and stromal scores, in line with the immune type and normal methylation type, demonstrating drug sensitivity to PD-1 inhibitors. This cluster was hardly engaged in metabolic signatures. The C3 subtype had the highest enrichment scores for carcinogenic pathways and the lowest enrichment scores for immune infiltration. It exhibited lower enrichment for metabolic signatures than C1 but greater enrichment than C2 for metabolic signatures. Overall, this research was undertaken to explore the metabolic landscape of SKCM and separately detect three clusters with different characteristics.
The results indicated that C1 had the most differential metabolic pathways, with 21 associated metabolic pathways all upregulated in this subtype. Therefore, we defined C1 as the metabolically active subtype. Several recent studies have focused on building a connection between endogenous tumor metabolism and immunotherapy. For example, recent studies have indicated that an increase in glycolytic metabolism in SKCM is linked to resistance to adoptive T-cell and immune checkpoint blockade therapies (Cascone et al., 2018). The hypermetabolic activation of C1 induces, on the one hand, nutrient depletion and hypoxia in the tumor microenvironment, thus, establishing metabolic competition between tumor cells and infiltrating immune cells. On the other hand, the active metabolism of SKCM in the tumor microenvironment leads to toxic concentrations of certain metabolites, including elevated concentrations of adenosine, kynurenine, ornithine, and reactive oxygen species, all of which markedly inhibit the anti-tumor immune response (Leone and Powell, 2020). The taxonomic classification model of SKCM reported in 2015 identified three transcriptomic clusters based on mRNA transcriptional gene function, namely, “immune,” “keratin,” and “MITF-low” (2015). Of the three subtypes, the keratin subtype had the poorest prognosis among patients with localized melanoma metastases and featured high expression of keratin-coding genes and metabolic genes in organ development, which was consistent with C1 having the poorest prognosis and characterized by enrichment gene signatures of keratin-related metabolic procedures. In addition, the clinical signature analysis showed that most patients in C1 were in the advanced pathological stage, followed by C3.
Several metabolic processes were upregulated in subclass C1, including amino acid metabolism, carbohydrate metabolism, and lipid metabolism. Enrichment results of metabolic signatures indicate the possibility that C1 could benefit from targeted metabolism therapy. Metabolic therapies targeting certain metabolic processes provide alternatives for chemoresistant patients. For instance, metformin can exert an anti-melanoma effect (Li et al., 2018) and promote combination treatment efficacy with anti-PD-1 and anti-cytotoxic T-lymphocyte-associated protein-4 (CTLA4) agents (Afzal et al., 2018; Cha et al., 2018). Identifying the potential beneficiaries of metabolic therapies has always proven to be challenging (Goodpaster and Sparks, 2017). Meanwhile, some key processes, such as glycolytic metabolism, hexosamine biosynthesis pathway, and glutathione metabolic pathways, were enriched in C3, suggesting that patients with this (HBP) are a shunt pathway of glycolysis and a key metabolic juncture in cancer cells. Recent studies have demonstrated that treatment with small molecule drugs that target HBP lead to increased infiltration of CD8 + T cells, sensitizing pancreatic tumors to anti-PD1 therapy and causing tumor regression and prolonged survival (Sharma et al., 2020). This implies that we may achieve precise targeting of C3 subclasses by targeting HBP in SKCM. This subcategory survey provides new clues to forecast potential candidates for metabolic therapies, which needs further verification in a large clinical cohort.
Tumor immune infiltration results showed that C2 had elevated levels of all immune cell lineages, indicating that it was in a highly immune-activating state. We believe that it contributes to the best prognosis of C2 among the three subtypes. C2 matches the immune subtype and is accompanied by the highest associated immune score. The immune subtype exhibited high expression of immune cell subsets (T cells, B cells, and NK cells), immune signaling molecules, and immune checkpoint-related genes, and patients with regional metastatic SKCM in the immune subtype had a better prognosis than the other two subtypes due to their active host immune response (2015). Furthermore, tumor microenvironment-related assessments revealed that C2 had higher immune and stromal scores. Research on SKCM classification based on TMB found that those with high immune scores had a favorable prognosis, and those with low immune scores had a poor prognosis (Hu et al., 2020), and these findings were consistent with the favorable prognosis of the immune-high C2 subtype in our study. Immune checkpoint blocking antibodies, such as pembrolizumab, nivolumab, and ipilimumab, which target PD-1 or CTLA-4, have considerably transformed the therapeutic landscape of SKCM in recent years (Schvartsman et al., 2019). C2 exhibited a higher expression of immune checkpoint genes, especially LAG3, CD274 (also known as PD-L1), PDCD1, CD247, CTLA4, PDCD1LG2, TNFRSF9, TNFRSF4, and TLR9, demonstrating that it shows promising sensitivity against anti-PD-1 therapy and other checkpoint inhibitors. The high expression enrichment score of PD-1 (CD274) may contribute to the susceptibility of C2 to anti-PD-1 therapy. This aligns with the results of the TIDE algorithm, in which C2 was sensitive to PD-1 immunotherapy. In contrast, immune checkpoint expression of fibrinogen-like protein 1 (FGL1) was highest in the C3 subtype, and FGL1 is the major LAG-3 functional ligand that acts independently of MHC-II (Wang et al., 2019). FGL1 represses antigen-specific T-cell activation and exerts a tumor immunosuppressive effect, and it is associated with poor prognosis and resistance to anti-PD-1 therapy.
The C3 subtype paralleled the MITF-low subtype of SKCM. The “MITF-low” cluster is marked with a low expression of genes associated with immune regulation (2015) and pigmentation markers, which is consistent with the low immune cell infiltration in C3. Meanwhile, studies have consistently reported that low MITF expression is an early resistance to multiple targeted drugs (Müller et al., 2014). Furthermore, the enrichment scores of C3 were significantly higher than those of C1 and C2 in most carcinogenic pathways, such as NRF2, PI3K, TGF-β, TP53, and Hippo, which accounts for the poor prognosis of patients in the C3 subgroup. In particular, activation of NRF2 inhibits the activity of the melanocyte lineage marker MITF and blunts the induction of innate immune responses in SKCM (Jessen et al., 2020). This is in line with the similarity of C3 with the MITF-low cluster. These data suggested that the C3 subclass was of high heterogeneity and might be refractory. Notably, previous studies demonstrated that the MITF-low subpopulation can be reversed by combining NK-κB inhibitors with SKCM resistance to BRAF inhibitors (Konieczkowski et al., 2014; Su et al., 2020), implying that the C3 subclass would best respond to the combination of NFK-κB inhibitors and BRAF inhibitors. C3 was also similar to the CpG island methylator phenotype (CIMP) (2015), one of the reported molecular subtypes of SKCM. Studies have shown that the CIMP pattern is implicated in the progression of the clinical stage of malignant melanoma (Tanemura et al., 2009). This is due to the fact that patients with SKCM of the C3 subtype are mostly at a higher pathological grade. Therefore, the poor prognosis in C3 may be attributed to the combined effects of low infiltration of immune cells, high enrichment of carcinogenesis pathways, low MITF expression, CIMP, and advanced pathological stage.
Skin cutaneous melanoma has an extremely high TMB due to ultraviolet (UV) mutagenesis (Lo et al., 2021), and these neoantigen burdens may alter T-cell responses in the tumor microenvironment (Hollern et al., 2019). Therefore, we examined whether copy number aberrations (deletions and amplifications), TMB, and neoantigens could be associated with the subtypes of SKCM. The results illustrated no statistically significant correlation between neoantigen burden and subtypes, with the TMB in C2 being the highest in the three types. Regarding copy number aberrations, C2 patients had fewer amplifications and deletions than C1 and C3. Our data revealed a high TMB in the C2 subclass, which is compatible with the result that the C2 subclass was sensitive to PD-1 treatment. Moreover, the number of somatic mutations was higher in C1 and C3 than in C2, which may lead to a poorer prognosis for C1 and C3. Notably, the high frequency of BRAF mutations in C2 and C3 implies that these two subtypes of patients may benefit from BRAF-targeted inhibitors (Menzer et al., 2019). In t-SNE analysis, C1 separated into two subpopulations, and some outliers of C1 were mixed with C3. We speculate that C1, which shows heterogeneity, is likely to be subdivided into two different subtypes, one of which is similar to C3 accompanied by similar stromal scores and mutation plots. In addition, C3 and C1 showed similar molecular patterns. Relatively higher metabolic signature level and worse prognosis, as well as lower abundance in immune signatures, distinguished C1 from C2. A more reasonable outcome will be attained if these C1 outliers are relabeled as C3.
Furthermore, we identified significant somatic mutation alteration sites in different subtypes according to GISTIC2.0. The number of loci and the degree of copy number variation of these genes were lower in C2 than in C1 and C3. In contrast, the number of amplification and deletion sites of the SKCM gene was higher in C3 than in C1 and C2. Furthermore, the magnitude of copy number variation showed the same trend. These results also showed a better prognosis of C2 than C1 and C3 and illustrated the association between the immune stromal infiltration of C1 and C3 and the occurrence of their tumor copy number variations. From the drug sensitivity analysis, we hypothesized that C2 may be more sensitive to chemotherapeutic agents than C1 and C3, suggesting that C2 may benefit from the combination of chemotherapy and immunotherapy. The GDSC drug sensitivity results indicated that two compounds, CGP-60474, a potent inhibitor of cyclin-dependent kinase (CDK), and bryostatin-1, a powerful protein kinase C (PKC) agonist, showed high drug sensitivity toward the three subtypes, representing promising chemotherapeutic drug candidates that target SKCM. The above discussion reveals that patients with the C2 subtype of SKCM may benefit from a combination of chemotherapy and immunotherapy. This study aimed to substantiate the need for personalized and precise treatment in clinical practice. The above discussion reveals that the three subtypes of SKCM may potentially benefit from the combination of chemotherapy and immunotherapy. This study aimed to provide supporting evidence for the availability of candidate target agents in clinical practice and achieve personalized precision therapy.
In the 90-gene classifier, many of the metabolic genes have been reported to play an important role in the inner metabolism and progression of SKCM, such as tyrosinase-related protein 1 (TYRP1) (Gilot et al., 2017) and OCA2 melanosomal transmembrane protein (OCN2) (Peterson et al., 2019), in particular, SLC7A4 of the solute carrier family of C1 to the transportation of cationic amino acids (Wolf et al., 2002). However, there have been no studies on SLC7A4 in SKCM. Therefore, we analyzed the expression signature of SLC7A4 and its prognostic value in TCGA-SKCM. We performed an IHC analysis of SLC7A4 from a clinical cohort to validate its role in SKCM. The results demonstrated that SLC7A4 was more highly expressed in melanoma tissues than in normal skin tissues, and the relatively higher expression of SLC7A4 presented a poor prognosis in SKCM patients, indicating that it could serve as a prognostic biomarker for SKCM. Moreover, the IHC results illustrated that the expression of SLC7A4 was proportional to the progression of melanoma. Therefore, SLC7A4 could be utilized as a therapeutic target for SKCM. These reports and experimental analyses reinforced the reliability of this gene classifier to a certain extent.
To the best of our knowledge, the present study is the first to generate a metabolic classification of SKCM that confirms the results of previous studies on SKCM subtypes while retaining their characteristics. Specifically, this classification matched the three subtypes from TCGA (keratin type, immune type, and MITF-low). C1 matched the keratin type and was highly metabolically active. C2 was the immune type, characterized by high immune infiltration, high expression of immune checkpoints, high immune and stromal scores, and favorable prognosis. We demonstrated for the first time that the C2 subtype corresponding to the immune type was prone to be a responder to PD-1 immunotherapy. C3 corresponded with MITF-low and was enriched in CpG island methylation. We also performed a preliminary experimental validation of the role of SLC7A4 in SKCM from the 90- gene classifier. Hence, the current study constitutes a novel demonstration of the presence of TCGA subtypes in TCGA-SKCM cohort and the GEO-SKCM cohort. Furthermore, this study not only validates the clinical significance of TCGA-SKCM classification but also reveals unexplored features, such as tumor microenvironment signatures and response to immunotherapy. This study provides new insights into the heterogeneity of SKCM in terms of the metabolic landscape by classifying SKCM into three clusters, active, intermediate, and depleted metabolic activity, which suggests possible therapeutic options for each subtype. However, we acknowledge several flaws in this study. First of all, more clinical data are needed to verify the reliability of this classification standard. Second, additional validation like cellular and molecular experiments of the classifier is required to prove the findings.
In conclusion, here, we proposed novel classifications of SKCM from a metabolic perspective with three subtypes, namely, metabolically active, intermediate, and depleted. C1 was closely associated with metabolic processes and had the worst prognosis, consistent with the keratin subtype. C2 exhibited higher immune infiltration, high immune and stromal scores, and sensitivity to PD-1 immune blockade agents, correlating with the immune type, and it had a favorable prognosis. C3 had an enriched carcinogenic pathway with a high degree of prognosis and relatively poor prognosis, and it was less metabolically active than C1 but more active than C2. Furthermore, the generated 90-gene classifier had a high predictive value for SKCM, and this classifier may help predict the prognosis of patients with SKCM more accurately and provide precise therapeutic approaches for these patients. Thus, our study further enhanced the recognition of the metabolic hallmarks of SKCM and contributed valuable new information regarding SKCM subtypes such that patients with this disease can receive personalized immunotherapy and more accurate prognosis prediction.
Data Availability Statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
Ethics Statement
The studies involving human participants were reviewed and approved by the Ethics Committee of Foshan First People’s Hospital. The patients/participants provided their written informed consent to participate in this study.
Author Contributions
RY and ZW conceived and designed the study. SZ, JL, ZW, RY, and YQ contributed to the data acquisition. JM, RG, and XP analyzed the data. ZW and SZ wrote the manuscript. All authors read and approved the final manuscript.
Funding
This study was supported by the National Natural Science Foundation of China (Nos. 82002913 and 81772136), the Foundation of Foshan City (Nos. FS0AA-KJ218-1301-0034 and 2018AB003411), the Special Fund of the Foshan Summit Plan (Nos. 2019C002, 2019D008, 2019A006, and 2020A015), the Foundation of Sichuan Science and Technology Department (2019YJ0652), Guang Dong Basic and Applied Basic Research Foundation (2021A1515011453), and the Medical Science Foundation of Cheng Du (2018003).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.707677/full#supplementary-material
Supplementary Figure 1 | GO enrichment analysis of differential genes by Goplot package. GO enrichment of differential genes between C1 clusters and C2 and C3 clusters by the Goplot package, which is visualized as a taxonomic bubble plot (A) and a circle plot (B). Venn diagram (C) of differential gene enrichment in “autophagy,” “protein binding,” “lysosome organization” three biological processes. Goplot package for GO enrichment of differential genes between C2 clusters and C1 and C3 clusters, and visualized as categorical bubble plots (D) and circle plots (E). Venn diagram (F) of differential gene enrichment in the three biological processes “immune response,” “regulation of immune response,” “adaptive immune response.” GO enrichment of differential genes between C3 clusters and C1 and C2 clusters was carried out by the Goplot package and visualized in the form of classification bubble plots (G) and circle plots (H). Venn diagram (I) of differential genes enriched in three biological processes: “protein binding,” “nucleotide-excision repair,” “DNA damage recognition,” and “poly(A) RNA binding.”
Supplementary Figure 2 | KEGG pathway enrichment analysis was performed using KOBAS 3.0. The pathway figures depicted the common KEGG enrichment pathways for genes that are different between the three clustered subtypes and other subtypes. p < 0.05 pathways are considered significantly enriched. Apoptosis (A), Cell cycle (B), HIF-1 signaling pathway (C), Human T-cell leukemia virus 1 (D), Pathways in cancer infection (E), PI3K-Akt signaling pathway (F), red color stands for enriched genes with up-regulated expression.
Supplementary Figure 3 | Mutation profiling of 90-gene classifier in the pan-cancer analysis and the TF-mRNA-miRNA network construction. (A) The gene mutation and gene copy number of 90-gene classifier in 32 TCGA pan-cancer databases were obtained from the cBioportal online database, and the analysis of the mutation was performed to illustrate the distribution proportion of each mutation type of this gene classifier in pan-cancer types, including gene mutation, fusion, amplification, deep deletion and multiple alterations in a bar chart. (B) Construction of potential transcription factor and miRNA (TF-mRNA-miRNA) interaction network for 90-gene classifier using Cytoscape software, with red representing transcription factors, blue representing mRNA, and green representing miRNA.
Supplementary Table 1 | The information of clinical characteristics of TCGA-SKCM cohorts.
Supplementary Table 2 | Independency test between meta cluster and other clinical information and previous classification in TCGA-SKCM cohorts.
Supplementary Table 3 | Independency test between meta cluster and other clinical information in GEO-SKCM cohorts.
Supplementary Table 4 | The ratio of driver-type oncogenic mutations based on TCGA-SKCM in each subtype.
Supplementary Table 5 | The gene list of the 90 gene-classifier.
Supplementary Table 6 | Analysis of GO enrichment of differential genes between subtypes C1 and C2 and C3 by Goplot package.
Supplementary Table 7 | GO enrichment analysis of differential genes between subtypes C2 and subtypes C1 and C3 by Goplot package.
Supplementary Table 8 | The GO enrichment analysis of differential genes between subtypes C3 and C1 and C2 by the Goplot package.
Supplementary Table 9 | KEGG pathway enrichment analysis was performed using KOBAS 3.0 (p < 0.05), and the table showed the Top10 enriched KEGG pathways.
Abbreviations
SKCM, skin cutaneous melanoma; ICI, immune checkpoint inhibitors; CTLA-4, cytotoxic T-lymphocyte-associated protein-4; GDSC, genomics of drug sensitivity in cancer; GEO, Gene Expression Omnibus; TCGA, The Cancer Genome Atlas; Tem, effective memory T cell; GSVA, gene set variation analysis; GSEA, gene set enrichment analysis; SLC7A4, solute carrier family 7 member 4; NMF, non-negative matrix factorization; CIMP, CpG island methylator phenotype; OS, overall survival.
Footnotes
- ^ http://firebrowse.org/
- ^ http://software.broadinstitute.org/gsea/msigdb/index.jsp
- ^ https://cran.r-project.org/web/packages/ggplot2/index.html
- ^ http://kobas.cbi.pku.edu.cn
- ^ https://xenabrowser.net/datapages/
- ^ http://www.genepattern.org/
References
Afzal, M. Z., Mercado, R. R., and Shirai, K. (2018). Efficacy of metformin in combination with immune checkpoint inhibitors (anti-PD-1/anti-CTLA-4) in metastatic malignant melanoma. J. Immunother. Cancer 6;64. doi: 10.1186/s40425-018-0375-1
Alexandrov, L. B., Kim, J., Haradhvala, N. J., Huang, M. N., Tian Ng, A. W., Wu, Y., et al. (2020). The repertoire of mutational signatures in human cancer. Nature 578, 94–101. doi: 10.1038/s41586-020-1943-3
Barbie, D. A., Tamayo, P., Boehm, J. S., Kim, S. Y., Moody, S. E., Dunn, I. F., et al. (2009). Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature 462, 108–112. doi: 10.1038/nature08460
Becht, E., Giraldo, N. A., Lacroix, L., Buttard, B., Elarouci, N., Petitprez, F., et al. (2016). Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol. 17:218. doi: 10.1186/s13059-016-1070-5
Bolick, N. L., and Geller, A. C. (2021). Epidemiology of melanoma. Hematol. Oncol. Clin. North Am. 35, 57–72. doi: 10.1016/j.hoc.2020.08.011
Bristot, I. J., Kehl Dias, C., Chapola, H., Parsons, R. B., and Klamt, F. (2020). Metabolic rewiring in melanoma drug-resistant cells. Crit. Rev. Oncol. Hematol. 153:102995. doi: 10.1016/j.critrevonc.2020.102995
Brunet, J. P., Tamayo, P., Golub, T. R., and Mesirov, J. P. (2004). Metagenes and molecular pattern discovery using matrix factorization. Proc. Natl. Acad. Sci. U.S.A. 101, 4164–4169. doi: 10.1073/pnas.0308531101
Cancer Genome Atlas Network (2015). Genomic classification of cutaneous melanoma. Cell 161, 1681–1696. doi: 10.1016/j.cell.2015.05.044
Cascone, T., McKenzie, J. A., Mbofung, R. M., Punt, S., Wang, Z., Xu, C., et al. (2018). Increased tumor glycolysis characterizes immune resistance to adoptive T Cell therapy. Cell Metab. 27, 977–987.e4. doi: 10.1016/j.cmet.2018.02.024
Cha, J. H., Yang, W. H., Xia, W., Wei, Y., Chan, L. C., Lim, S. O., et al. (2018). Metformin promotes antitumor immunity via endoplasmic-reticulum-associated degradation of PD-L1. Mol. Cell 71, 606–620.e7. doi: 10.1016/j.molcel.2018.07.030
Colaprico, A., Silva, T. C., Olsen, C., Garofano, L., Cava, C., Garolini, D., et al. (2016). TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res. 44:e71. doi: 10.1093/nar/gkv1507
Gao, J., Aksoy, B. A., Dogrusoz, U., Dresdner, G., Gross, B., Sumer, S. O., et al. (2013). Integrative analysis of complex cancer genomics and clinical profiles using the cBioPortal. Sci. Signal. 6:l1. doi: 10.1126/scisignal.2004088
Gaujoux, R., and Seoighe, C. (2010). A flexible R package for nonnegative matrix factorization. BMC Bioinformatics 11:367. doi: 10.1186/1471-2105-11-367
Ghiorzo, P., Gargiulo, S., Pastorino, L., Nasti, S., Cusano, R., Bruno, W., et al. (2006). Impact of E27X, a novel CDKN2A germ line mutation, on p16 and p14ARF expression in Italian melanoma families displaying pancreatic cancer and neuroblastoma. Hum. Mol. Genet. 15, 2682–2689. doi: 10.1093/hmg/ddl199
Gibcus, J. H., Kok, K., Menkema, L., Hermsen, M. A., Mastik, M., Kluin, P. M., et al. (2007). High-resolution mapping identifies a commonly amplified 11q13.3 region containing multiple genes flanked by segmental duplications. Hum. Genet. 121, 187–201. doi: 10.1007/s00439-006-0299-6
Gilot, D., Migault, M., Bachelot, L., Journé, F., Rogiers, A., Donnou-Fournet, E., et al. (2017). A non-coding function of TYRP1 mRNA promotes melanoma growth. Nat. Cell Biol. 19, 1348–1357. doi: 10.1038/ncb3623
Goodpaster, B. H., and Sparks, L. M. (2017). Metabolic flexibility in health and disease. Cell Metab. 25, 1027–1036. doi: 10.1016/j.cmet.2017.04.015
Guy, G. P., Thomas, C. C., Thompson, T., Watson, M., Massetti, G. M., and Richardson, L. C. (2015). Vital signs: melanoma incidence and mortality trends and projections - United States, 1982-2030. MMWR Morb. Mortal Wkly. Rep. 64, 591–596.
Hänzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 14:7. doi: 10.1186/1471-2105-14-7
Harel, M., Ortenberg, R., Varanasi, S. K., Mangalhara, K. C., Mardamshina, M., Markovits, E., et al. (2019). Proteomics of melanoma response to immunotherapy reveals mitochondrial dependence. Cell 179, 236–250.e18. doi: 10.1016/j.cell.2019.08.012
Havel, J. J., Chowell, D., and Chan, T. A. (2019). The evolving landscape of biomarkers for checkpoint inhibitor immunotherapy. Nat. Rev. Cancer 19, 133–150. doi: 10.1038/s41568-019-0116-x
Hollern, D. P., Xu, N., Thennavan, A., Glodowski, C., Garcia-Recio, S., Mott, K. R., et al. (2019). B cells and T follicular helper cells mediate response to checkpoint inhibitors in high mutation burden mouse models of breast cancer. Cell 179, 1191–1206.e21. doi: 10.1016/j.cell.2019.10.028
Hoshida, Y., Brunet, J. P., Tamayo, P., Golub, T. R., and Mesirov, J. P. (2007). Subclass mapping: identifying common subtypes in independent disease data sets. PLoS One 2:e1195. doi: 10.1371/journal.pone.0001195
Hu, B., Wei, Q., Li, X., Ju, M., Wang, L., Zhou, C., et al. (2020). Development of an IFNγ response-related signature for predicting the survival of cutaneous melanoma. Cancer Med. 9, 8186–8201. doi: 10.1002/cam4.3438
Huang da, W., Sherman, B. T., and Lempicki, R. A. (2009). Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. 37, 1–13. doi: 10.1093/nar/gkn923
Huang, D. W., Sherman, B. T., Tan, Q., Kir, J., Liu, D., Bryant, D., et al. (2007). DAVID bioinformatics resources: expanded annotation database and novel algorithms to better extract biology from large gene lists. Nucleic Acids Res. 35, W169–W175. doi: 10.1093/nar/gkm415
Jessen, C., Kreß, J. K. C., Baluapuri, A., Hufnagel, A., Schmitz, W., Kneitz, S., et al. (2020). The transcription factor NRF2 enhances melanoma malignancy by blocking differentiation and inducing COX2 expression. Oncogene 39, 6841–6855. doi: 10.1038/s41388-020-01477-8
Konieczkowski, D. J., Johannessen, C. M., Abudayyeh, O., Kim, J. W., Cooper, Z. A., Piris, A., et al. (2014). A melanoma cell state distinction influences sensitivity to MAPK pathway inhibitors. Cancer Discov. 4, 816–827. doi: 10.1158/2159-8290.Cd-13-0424
Leone, R. D., and Powell, J. D. (2020). Metabolism of immune cells in cancer. Nat. Rev. Cancer 20, 516–531. doi: 10.1038/s41568-020-0273-y
Li, K., Zhang, T. T., Wang, F., Cui, B., Zhao, C. X., Yu, J. J., et al. (2018). Metformin suppresses melanoma progression by inhibiting KAT5-mediated SMAD3 acetylation, transcriptional activity and TRIB3 expression. Oncogene 37, 2967–2981. doi: 10.1038/s41388-018-0172-9
Lo, J. A., Kawakubo, M., Juneja, V. R., Su, M. Y., Erlich, T. H., LaFleur, M. W., et al. (2021). Epitope spreading toward wild-type melanocyte-lineage antigens rescues suboptimal immune checkpoint blockade responses. Sci. Transl. Med. 13:eabd8636. doi: 10.1126/scitranslmed.abd8636
Lu, X., Meng, J., Zhou, Y., Jiang, L., and Yan, F. (2020). MOVICS: an R package for multi-omics integration and visualization in cancer subtyping. Bioinformatics 36, 5539–5541. doi: 10.1093/bioinformatics/btaa1018
Mayakonda, A., Lin, D. C., Assenov, Y., Plass, C., and Koeffler, H. P. (2018). Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 28, 1747–1756. doi: 10.1101/gr.239244.118
Menzer, C., Menzies, A. M., Carlino, M. S., Reijers, I., Groen, E. J., Eigentler, T., et al. (2019). Targeted therapy in advanced melanoma with rare BRAF mutations. J. Clin. Oncol. 37, 3142–3151. doi: 10.1200/jco.19.00489
Mermel, C. H., Schumacher, S. E., Hill, B., Meyerson, M. L., Beroukhim, R., and Getz, G. (2011). GISTIC2.0 facilitates sensitive and confident localization of the targets of focal somatic copy-number alteration in human cancers. Genome Biol. 12:R41. doi: 10.1186/gb-2011-12-4-r41
Müller, J., Krijgsman, O., Tsoi, J., Robert, L., Hugo, W., Song, C., et al. (2014). Low MITF/AXL ratio predicts early resistance to multiple targeted drugs in melanoma. Nat. Commun. 5:5712. doi: 10.1038/ncomms6712
Nassar, K. W., and Tan, A. C. (2020). The mutational landscape of mucosal melanoma. Semin. Cancer Biol. 61, 139–148. doi: 10.1016/j.semcancer.2019.09.013
Pavlova, N. N., and Thompson, C. B. (2016). The emerging hallmarks of cancer metabolism. Cell Metab. 23, 27–47. doi: 10.1016/j.cmet.2015.12.006
Peifer, M., Hertwig, F., Roels, F., Dreidax, D., Gartlgruber, M., Menon, R., et al. (2015). Telomerase activation by genomic rearrangements in high-risk neuroblastoma. Nature 526, 700–704. doi: 10.1038/nature14980
Peterson, K. A., Neuffer, S., Bean, M. E., New, L., Coffin, A. B., and Cooper, C. D. (2019). Melanosome maturation proteins Oca2, Mitfa and Vps11 are differentially required for cisplatin resistance in zebrafish melanocytes. Exp. Dermatol. 28, 795–800. doi: 10.1111/exd.13937
Possemato, R., Marks, K. M., Shaul, Y. D., Pacold, M. E., Kim, D., Birsoy, K., et al. (2011). Functional genomics reveal that the serine synthesis pathway is essential in breast cancer. Nature 476, 346–350. doi: 10.1038/nature10350
Ratnikov, B. I., Scott, D. A., Osterman, A. L., Smith, J. W., and Ronai, Z. A. (2017). Metabolic rewiring in melanoma. Oncogene 36, 147–157. doi: 10.1038/onc.2016.198
Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 43:e47. doi: 10.1093/nar/gkv007
Robert, C., Thomas, L., Bondarenko, I., O’Day, S., Weber, J., Garbe, C., et al. (2011). Ipilimumab plus dacarbazine for previously untreated metastatic melanoma. N. Engl. J. Med. 364, 2517–2526. doi: 10.1056/NEJMoa1104621
Roh, W., Chen, P. L., Reuben, A., Spencer, C. N., Prieto, P. A., Miller, J. P., et al. (2017). Integrated molecular analysis of tumor biopsies on sequential CTLA-4 and PD-1 blockade reveals markers of response and resistance. Sci. Transl. Med. 9:eaah3560. doi: 10.1126/scitranslmed.aah3560
Rosario, S. R., Long, M. D., Affronti, H. C., Rowsam, A. M., Eng, K. H., and Smiraglia, D. J. (2018). Pan-cancer analysis of transcriptional metabolic dysregulation using The Cancer Genome Atlas. Nat. Commun. 9:5330. doi: 10.1038/s41467-018-07232-8
Ruocco, M. R., Avagliano, A., Granato, G., Vigliar, E., Masone, S., Montagnani, S., et al. (2019). Metabolic flexibility in melanoma: a potential therapeutic target. Semin. Cancer Biol. 59, 187–207. doi: 10.1016/j.semcancer.2019.07.016
Sanchez-Vega, F., Mina, M., Armenia, J., Chatila, W. K., Luna, A., La, K. C., et al. (2018). Oncogenic signaling pathways in The Cancer Genome Atlas. Cell 173, 321–337.e10. doi: 10.1016/j.cell.2018.03.035
Schadendorf, D., Hodi, F. S., Robert, C., Weber, J. S., Margolin, K., Hamid, O., et al. (2015). Pooled analysis of long-term survival data from phase II and phase III trials of ipilimumab in unresectable or metastatic melanoma. J. Clin. Oncol. 33, 1889–1894. doi: 10.1200/jco.2014.56.2736
Schvartsman, G., Taranto, P., Glitza, I. C., Agarwala, S. S., Atkins, M. B., and Buzaid, A. C. (2019). Management of metastatic cutaneous melanoma: updates in clinical practice. Ther. Adv. Med. Oncol. 11:1758835919851663. doi: 10.1177/1758835919851663
Sharma, N. S., Gupta, V. K., Garrido, V. T., Hadad, R., Durden, B. C., Kesh, K., et al. (2020). Targeting tumor-intrinsic hexosamine biosynthesis sensitizes pancreatic cancer to anti-PD1 therapy. J. Clin. Invest. 130, 451–465. doi: 10.1172/jci127515
Su, Y., Ko, M. E., Cheng, H., Zhu, R., Xue, M., Wang, J., et al. (2020). Multi-omic single-cell snapshots reveal multiple independent trajectories to drug tolerance in a melanoma cell line. Nat. Commun. 11:2345. doi: 10.1038/s41467-020-15956-9
Tanemura, A., Terando, A. M., Sim, M. S., van Hoesel, A. Q., de Maat, M. F., Morton, D. L., et al. (2009). CpG island methylator phenotype predicts progression of malignant melanoma. Clin. Cancer Res. 15, 1801–1807. doi: 10.1158/1078-0432.Ccr-08-1361
Thanasai, J., Limpaiboon, T., Jearanaikoon, P., Bhudhisawasdi, V., Khuntikeo, N., Sripa, B., et al. (2006). Amplification of D22S283 as a favorable prognostic indicator in liver fluke related cholangiocarcinoma. World J. Gastroenterol. 12, 4338–4344. doi: 10.3748/wjg.v12.i27.4338
Vivian, J., Rao, A. A., Nothaft, F. A., Ketchum, C., Armstrong, J., Novak, A., et al. (2017). Toil enables reproducible, open source, big biomedical data analyses. Nat. Biotechnol. 35, 314–316. doi: 10.1038/nbt.3772
Walter, W., Sánchez-Cabo, F., and Ricote, M. (2015). GOplot: an R package for visually combining expression data with functional analysis. Bioinformatics 31, 2912–2914. doi: 10.1093/bioinformatics/btv300
Wang, J., Sanmamed, M. F., Datar, I., Su, T. T., Ji, L., Sun, J., et al. (2019). Fibrinogen-like protein 1 is a major immune inhibitory ligand of LAG-3. Cell 176, 334–347.e12. doi: 10.1016/j.cell.2018.11.010
Wolf, S., Janzen, A., Vékony, N., Martiné, U., Strand, D., and Closs, E. I. (2002). Expression of solute carrier 7A4 (SLC7A4) in the plasma membrane is not sufficient to mediate amino acid transport activity. Biochem. J. 364(Pt 3), 767–775. doi: 10.1042/bj20020084
Woods, K., Knights, A. J., Anaka, M., Schittenhelm, R. B., Purcell, A. W., Behren, A., et al. (2016). Mismatch in epitope specificities between IFNγ inflamed and uninflamed conditions leads to escape from T lymphocyte killing in melanoma. J. Immunother. Cancer 4:10. doi: 10.1186/s40425-016-0111-7
Wu, J., Mao, X., Cai, T., Luo, J., and Wei, L. (2006). KOBAS server: a web-based platform for automated annotation and pathway identification. Nucleic Acids Res. 34, W720–W724. doi: 10.1093/nar/gkl167
Yang, W., Soares, J., Greninger, P., Edelman, E. J., Lightfoot, H., Forbes, S., et al. (2013). Genomics of drug sensitivity in cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells. Nucleic Acids Res. 41, D955–D961. doi: 10.1093/nar/gks1111
Yoshihara, K., Shahmoradgoli, M., Martínez, E., Vegesna, R., Kim, H., Torres-Garcia, W., et al. (2013). Inferring tumour purity and stromal and immune cell admixture from expression data. Nat. Commun. 4:2612. doi: 10.1038/ncomms3612
Yu, G., Wang, L. G., Han, Y., and He, Q. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. Omics 16, 284–287. doi: 10.1089/omi.2011.0118
Zhang, Y., Parmigiani, G., and Johnson, W. E. (2020). ComBat-seq: batch effect adjustment for RNA-seq count data. NAR Genom. Bioinform. 2:lqaa078. doi: 10.1093/nargab/lqaa078
Zhou, G., Soufan, O., Ewald, J., Hancock, R. E. W., Basu, N., and Xia, J. (2019). NetworkAnalyst 3.0: a visual analytics platform for comprehensive gene expression profiling and meta-analysis. Nucleic Acids Res. 47, W234–W241. doi: 10.1093/nar/gkz240
Keywords: metabolism subtypes, tumor microenvironment, skin cutaneous melanoma, immune signature, mutation landscape, immunotherapy response
Citation: Yang R, Wang Z, Li J, Pi X, Gao R, Ma J, Qing Y and Zhou S (2021) The Identification of the Metabolism Subtypes of Skin Cutaneous Melanoma Associated With the Tumor Microenvironment and the Immunotherapy. Front. Cell Dev. Biol. 9:707677. doi: 10.3389/fcell.2021.707677
Received: 10 May 2021; Accepted: 06 July 2021;
Published: 12 August 2021.
Edited by:
Na Luo, Nankai University, ChinaReviewed by:
Bingjie Li, First Affiliated Hospital of Zhengzhou University, ChinaYan Li, Nanjing Medical University, China
Copyright © 2021 Yang, Wang, Li, Pi, Gao, Ma, Qing and Zhou. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Yi Qing, qywxbb@hotmail.com; Sitong Zhou, sitongzhou@hotmail.com
†These authors have contributed equally to this work