Crosstalk of four kinds of cell deaths defines subtypes of cutaneous melanoma for precise immunotherapy and chemotherapy

Background: Cell death patterns can give therapeutic and biological clues that facilitate the development of individualized treatments for this lethal form of skin cancer.

Methods: We employed unsupervised clustering to establish robust classifications based on the four kinds of cell death-associated gene expression of 462 melanoma patients in the Cancer Genome Atlas (TCGA) and tested their reproducibility in two independent melanoma cohorts of 558 patients. We then used dimensionality reduction of graph learning to display the different characteristics of cell death patterns and immune microenvironments.

Results: We examined 570 cell death-associated gene expression data of melanoma patients for exploration, independent verification, and comprehensive classification of five reproducible melanoma subtypes (CS1 to CS5) with different genomic and clinical features. Patients in death-inactive subtypes (CS1, CS2, and CS5) had the least immune and stromal cell infiltration, and their prognosis was the poorest. A death-active subtype (CS4), on the other hand, had the highest infiltrated immune and stromal cells and elevated immune-checkpoints. As a result, these patients had the highest response to immunotherapy and the best prognosis. An additional subtype (CS3) had more diversified cell death and immune characteristics with moderate prognoses. Based on graph learning, we successfully divided the CS3 subtype into two subgroups (group A and group B) with distinct survival outcomes and immune features. Finally, we identified eight potential chemical drugs that were specifically targeted for the therapy of melanoma subtypes.

Conclusions: This research defines the intrinsic subtypes of melanoma based on the crosstalk of four kinds of cell deaths, which affords a blueprint for clinical strategies and guiding precise immunotherapy and chemotherapy for melanoma patients.

1. Introduction

Cutaneous melanoma is the most lethal type of skin cancer and is caused by the malignant alteration of melanocytes (1). The prevalence of melanoma has steadily risen during the last few decades (2). Therapy procedures have considerably improved in recent years as a result of a better knowledge of the major oncogenes and signaling pathways involved in its development and progression—for instance, molecularly targeted therapies such as the inhibitors of BRAF and MEK and immune-checkpoint blockade like anti-PD1 and CTLA-4 treatments have now significantly improved survival in these melanoma patients (3, 4). However, responses to these therapies are still inconsistent. Many individuals do not benefit at all or relapse after a brief time of remission. Therefore, future results in these melanomas will be dependent on the discovery of new treatment targets and strategies to improve the present targets and immunotherapies.

Cell death is significant in the formation of the organism and sustains homeostasis to avoid disease development (5). Melanoma growth, like other types of cancer, is influenced by several cell death mechanisms and tumor microenvironments (6). The diverse microenvironments of cell death processes influence the immune response owing to the changes in tumor cell death and immune, and stromal cell activities. Thus, treatment efforts should focus on various cell death processes. Cell death is traditionally classified into programmed and non-programmed types depending on the modulation of the underlying processes (7). There are three classical types of programmed cell death (PCD) including apoptosis, necroptosis, and pyroptosis. Apoptosis belongs to non-lytic cell death, while necroptosis and pyroptosis mainly refer to lytic cell death (8, 9). These types of cell death cause intracellular component leakage, containing damage-associated molecular pattern molecules, which trigger a significant inflammatory response, also defined as inflammatory death (10). Furthermore, ferroptosis is a novel cell death based on iron regulation, which belongs to PCD but has unique biochemical and morphological variations from other PCDs such as pyroptosis and necroptosis (11, 12). Cell death has been found to influence the development of a variety of chronic illnesses—for example, dysfunction of the apoptosis pathway has a significant role in resistance to traditional anticancer medicines like targeted therapy, chemotherapy, and radiotherapy (13). Ferroptosis dysregulation is connected to carcinogenesis and has been proven to link the formation of melanoma, breast, gastric, and lymphoma (14). In addition, several research published in the last 5 years have revealed that tumor cells undergoing necroptosis, pyroptosis, and ferroptosis can induce a powerful anti-cancer immunity in in vitro and in vivo conditions, and their efficacy can be cohesively enhanced by checkpoint blockades, even in cancers with immunotherapy resistance (15–17).

Melanoma is one of the most diverse cancers, with varying degrees of aggressiveness among subtypes that necessitate different treatment strategies. Although stimulation of a specific form of cell death has progressively developed as a new strategy treatment for cancers, most cancers have a built-in resistance to a certain kind of cell death. Therefore, the crosstalk of four kinds of cell deaths (“apoptosis”, “ferroptosis”, “necroptosis”, and “pyroptosis”) has been discussed systematically in this research. Based on cell death-related gene signatures, melanoma patients were successfully classified into distinct subtypes with different clinical characteristics, molecular features, and responses to immunotherapy and chemotherapy.

2. Materials and methods

2.1. Melanoma collection and process

Publicly accessed databases like Gene-Expression Omnibus (GEO) and the Cancer Genome Atlas (TCGA) were used to retrieve the RNA-seq of melanoma data as well as corresponding clinical annotations. Following data processing steps such as converting to TPMs format, removing low expression values or missing values, and log2-transformation, three melanoma cohorts including the cutaneous melanoma of TCGA data (TCGA_SKCM), meta-GEO data, and meta-immune response data were selected for subsequent analysis. The TCGA_SKCM cohort was used as a training dataset for melanoma classification. Meanwhile, meta-GEO and meta-immune response data were regarded as testing datasets to prove the accuracy of classification. The meta-GEO cohort consisted of four GEO melanoma datasets including GSE19234, GSE53118, GSE54467, and GSE65904. Moreover, the meta-immune response cohort consists of the respective study of Hugo et al. (GSE78220) (18), Riaz et al. (GSE91061) (19), and Van Allen et al. (20), which has detailed information about clinical immunotherapy response. The “ComBat” algorithm of the sva Package was applied to merge several individual studies into the large meta-cohort. This method can reduce batch effects produced by distinct platforms and technological biases.

2.2. Four kinds of cell death-related gene

The Molecular Signatures Database in the Gene Set Enrichment Analysis (GSEA) website (https://www.gsea-msigdb.org) and previous literature in the PubMed database were used to compile lists of cell death-related genes. The search keywords include “apoptosis”, “ferroptosis”, “necroptosis”, “pyroptosis”, and “Homo sapiens”. After removing the overlapped genes, the list gene sets eventually included the apoptosis-related gene set (n = 247 unique genes), ferroptosis-related gene set (n = 206 unique genes), necroptosis-related gene set (n = 65 unique genes), and pyroptosis-related gene set (n = 52 unique genes) (Supplementary Table S1).

2.3. Clustering of four cell death expression patterns

In order to conduct systematic clustering, the transcriptome data of TCGA_SKCM were obtained to construct four cell death-related gene matrixes with rows corresponding to features and columns corresponding to common samples. Next, unsupervised clustering was applied to discover the potentially significant subtype of melanoma. We employed the “MOVICS” package to perform integrative classification and illustration for cancer subtyping research. The “MOVICS” package afforded two main cluster algorithms (gap statistics and clustering prediction index) to estimate the most possible number of clusters in a partition clustering (21). Therefore, we calculated the gap statistics and clustering prediction index (CPI) to find the best clustering number. As a result, 10 state-of-the-art integrative clustering approaches (CIMLR, iClusterBayes, COCA, ConsensusClustering, MoCluster, NEMO, IntNMF, LRAcluster, PINSPlus, and SNF) were independently performed to systematically cluster the TCGA_SKCM cohort. To boost clustering resilience, we adopted the notion of a consensus ensemble for subsequently classified results produced from diverse methods.

2.4. Tumor microenvironment calculation

Two methods for decoding microenvironment cells (ESTIMATE and Stemness Index Workflow) were adjusted to create a compendium of gene lists connected to particular microenvironment cells. The “ESTIMATE” approach was applied to estimate the infiltrated immune cells, stromal cells, and tumor purity in tumor tissue. The Stemness Index Workflow (https://bioinformaticsfmrp.github.io/PanCanStem_Web/) applied the one-class algorithm to calculate tumor stem cells in the tumor sample. We first employed these two bioinformatic methods to evaluate tumor immune indices (immune score, stromal score, tumor purity, and ESTIMATE score) and stem cell index (mRNAsi score) and then standardized the value of indices ranging from 0 to 1.

2.5. Exploring the distinct clinical, molecular, and cellular features of subtypes

Based on clinical features such as age, stage, metastasis, Clark level, and so on, we first assessed the different distribution of clinical features among the cell death subtypes in TCGA_SKCM cohort. To be specific, melanoma patients may benefit from subtype-specific mutations as a treatment target. Therefore, we compared the different mutational frequencies, tumor mutation burden (TMB), and fraction genome altered (FGA) among subtypes of cancer. Moreover, to identify subtype-specific functional pathways, the differentially expressed genes in each subtype were annotated in terms of Gene Ontology (GO) biological processes (c5.bp.v7.1.symbols.gmt) by Gene Set Enrichment Analysis (GSEA). Furthermore, we referred to 10 tumor-associated gene sets, five immune-microenvironment gene signatures, and seven kinds of immune-checkpoint molecules which were collected and deposited in the “IOBR” package. Then, we used the single-sample gene set enrichment analysis (ssGSEA) method to estimate the enrichment scores among melanoma subtypes.

2.6. Chemotherapy response prediction

The Genomics of Drug Sensitivity in Cancer (GDSC; https://www.cancerrxgene.org/) database was used to explore the unique genomics of drug sensitivity across melanoma cell subtypes for identifying candidate agents that displayed variable effectiveness in cell death-related gene categorized clusters. In this database, over 1,000 genetically defined human cancer cell lines have been treated with a range of anti-cancer drugs (367 compounds). In total, 34 kinds of these cell lines are derived from cutaneous melanoma. We first used KNN machine learning to categorize these cell lines into corresponding cell death gene-classified subtypes and then examined the variations in the area under the curve (AUC) value of drug responses among these subtypes. After removing drugs with more than 20% missing values, only 219 compounds were selected for further analysis. To evaluate if there is a significant difference, Kruskal–Wallis test with a p-value of 0.05 was utilized.

2.7. Effect of drugs on A357 melanoma cell

2.7.1. Chemicals

All chemicals were purchased from MedChemExpress (MCE, China) including ACY-1215 (no. HY-16026), tubastatin A (no. HY-13271A), and EHT 1864 (no. HY-16659).

2.7.2. Cell lines and cell culture

The human melanoma cell line A375 was a gift from Associate Professor Naihong Yan from the Research Laboratory of Ophthalmology, West China Hospital, Sichuan University, China. All cells were cultured in Dulbecco’s modified Eagle’s medium (Gibco, USA), supplemented with 10% fetal bovine serum (Gibco, USA) and 1% penicillin/streptomycin mixture. The cells were maintained in a humidified incubator (37°C, 5% CO₂, and 95% air).

2.7.3. Cell viability assay

Cell viability was detected with cell counting kit-8 (CCK8, MCE, China). A375 cells were plated in 96-well plates at a density of 1 × 104 cells per well and treated with different concentrations of ACY-1215, tubastatin A, and EHT 1864 (total volume of 200 μl per well). After 24 and 48 h, 10 μl of CCK8 solution was added to each well of the plate and then incubated at 37°C in the dark for 1 h. The optical density was measured at 450 nm by a microplate reader (Bio-Rad, USA).

2.7.4. Wound healing assay

The cells were seeded on six-well plates. After the cells have reached 100% confluence, the monolayer was scratched using a tip and washed with serum-free medium to remove detached cells. IC50 concentrations of ACY-1215, tubastatin A, and EHT 1864 were added to a serum-free medium for culture and photographed at 0, 24, and 48 h, respectively. The image-j software was applied to calculate the area of the wound.

2.7.5. Cell apoptosis assay

Annexin V-FITC/propidium iodide (PI) Apoptosis kit (BD Biosciences, USA) was used to measure the apoptosis level of A375 cells treated with IC50 concentration of ACY-1215, tubastatin A, and EHT 1864, respectively. The cells were stained with Annexin V and PI according to the instructions and then analyzed by flow cytometry (Beckman, USA). Early apoptotic A375 cells showed annexin V-positive and PI-negative, and late apoptotic A375 cells were stained with both annexin V- and PI-positive. The results were presented as the percentage of apoptotic cells.

2.8. Pseudotime analysis of subtypes

The pseudotime analysis based on dimensionality reduction of graph learning was carried out by applying the reduced dimension function of the “Monocle” package. The dimension reduction approach was discriminative dimensionality reduction with trees, and we set the maximum number of components at five. Finally, the melanoma patients were landscaped and plotted using the plot cell trajectory function, with the color correlating to the cell death subtypes mentioned above.

2.9. Statistical analysis

Statistical analysis was entirely carried out by using the R package (v.4.0.3) and related packages. Unsupervised clustering was conducted by “MOVICS” package (21). The “survival” and “survivalROC” packages were used for survival analysis. ESTIMATE and molecular signal pathways were estimated and visualized by “IOBR” package (22). Pseudotime analysis was visualized by the “Monocle” package. The Spearman test was used to assess the correlation coefficient. For comparisons of two or more groups, the Wilcoxon and Kruskal–Wallis tests were used. Fisher’s exact or chi-square test was used to examine the relationship between subgroup and clinicopathological features.

3. Results

3.1. Melanoma collection and clustering

After removing duplicated samples or samples without survival information, 462 melanoma patients were included in the TCGA_SKCM cohort. The meta-GEO cohort consisted of 412 melanoma patients, and the meta-immune response cohort consisted of 146 melanoma patients who received immune-checkpoint blockade therapy. A total of 570 genes that referred to four kinds of cell deaths (“apoptosis”, “ferroptosis”, “necroptosis”, and “pyroptosis”) were evaluated to explore the important clusters in the TCGA_SKCM cohort. In order to prove that these genes are consistent with cell death, we conducted GO enrichment analysis and found that these genes are positively activated in cell death-associated signal pathways such as regulation of the cysteine-type endopeptidase activity involved in the apoptotic process, regulation of apoptotic signaling pathway, and execution phase of apoptosis (Supplementary Figure S1). Next, CPI and gaps statistics uncovered that the optimum number of clusters k for clustering was three to five, and combined with the popular Clark level classifier for cutaneous melanoma, it has five classifications (Figure 1A). It should be highlighted that both CPI and gap statistics do not fall too much at k = 5. As a result of these considerations, k of 5 is selected as the ideal clustering number for further analysis. A consensus ensemble derived from 10 different clustering methods identified five robust subtypes of melanoma (Figure 1B). We defined these five clusters from CS1 to CS5. The complex cross-talk of four kinds of cell death expression may have a particular biological relevance that contributes to the heterogeneity of melanomas. Therefore, we generated a comprehensive heat map to present distinctive molecular patterns across apoptosis, ferroptosis, necroptosis, and pyroptosis (Figure 1C).

FIGURE 1

Figure 1 Identification of melanoma subtypes based on cell death-associated gene expression. (A) Calculating clustering prediction index (blue line) and Gaps-statistics (red line) in the cutaneous melanoma of The Cancer Genome Atlas data (TCGA_SKCM) cohort to determine the best cluster number. (B) Consensus heat map according to the findings of 10 integrative clustering methods with a cluster number of five. (C) Comprehensive heat map of distinctive molecular patterns across apoptosis, ferroptosis, necroptosis, and pyroptosis with annotation of genes. (D) Top 20 mutant genes in the TCGA_SKCM cohort. (E) Mutational OncoPrint of five identified melanoma subtypes in the TCGA_SKCM cohort.

3.2. Somatic mutation in subtypes

To identify the subtype-specific mutations in melanoma, the “maftools” package (23) was first employed to produce oncoPrint plots which illustrated the top 20 popular mutant genes including TTN (72%), MUC16 (67%), BRAF (50%), DNAH5 (49%), PCLO (45%) … (Figure 1D). The independent testing (Fisher’s exact or chi-square test) determined that seven genes (TTN, MUC16, BRAF, DNAH5, LRP1B, PKHD1L1, and APOB) have different mutant frequencies among these subtypes (Figure 1E). Among the frequently mutated genes, CS2 significantly harbored more mutations of TTN (86.5%; adj-P = 0.039), DNAH5 (54.1%; adj-P = 0.039), and APOB (37.8%; adj-P = 0.039) than the other subtypes, while CS3 and CS4 were significantly enriched in the mutations of BRAF (CS3: 57.4%; CS4: 60.2%; adj-P = 0.039), MUC16 (CS3: 76.1%; CS4: 71.1%; adj-P = 0.039), and LRP1B (48.3%; adj-P = 0.039) and were frequently mutated in CS3, whereas PKHD1L1 was significantly more mutated in CS2 than in other subtypes (40.5%; adj-P = 0.039) (Table 1). The FGA scores were then used to assess chromosomal instability, and we discovered that CS4 had superior chromosomal stability than the other subtypes, with a much lower fraction of genome lost or gained (Figure 2A). Moreover, TMB has arisen as an emerging biomarker for the prediction of numerous tumor types, prognosis, and immunotherapy responses. Therefore, TMB was contrasted among subtypes, and the Kruskal–Wallis test indicated that CS4 and CS3 have relatively higher TMB than the other subtypes (Figure 2B).

TABLE 1

Table 1 Independent test between subtype and mutation.

FIGURE 2

Figure 2 Independent verification of melanoma subtypes. (A) Bar plot of fraction genome altered among five identified melanoma subtypes in the cutaneous melanoma of The Cancer Genome Atlas data (TCGA_SKCM cohort). (B) Comparison of tumor mutation burden and single-nucleotide variants among five melanoma subtypes in the TCGA_SKCM cohort. (C) Kaplan–Meier survival curve of five melanoma subtypes in the TCGA_SKCM cohort. (D) Heat map of subtype-specific upregulated biomarkers for five melanoma subtypes in the TCGA_SKCM cohort. (E) Heat map of nearest template prediction in the meta-Gene Expression Omnibus (GEO) cohort using subtype-specific upregulated biomarkers identified from the TCGA_SKCM cohort. (F) Kaplan–Meier survival curve of five melanoma subtypes in the meta-GEO cohort. **** means p < 0.0001.

3.3. Survival analysis and validation of meta-GEO cohort

In the TCGA_SKCM cohort, 118 melanoma patients were classified into the CS1 subgroup, CS2 had 39 melanoma patients, and the CS3, CS4, and CS5 subgroups respectively consisted of 180, 84, and 42 patients. Among these subgroups, CS4 has the best survival prognosis than the other subgroups (overall survival, P < 0.001) (Figure 2C). The stratified survival analysis revealed that CS4 has a longer survival time than the patients in the CS1 (P = 0.002), CS2 (P < 0.001), and CS5 (P < 0.001) subtypes. Furthermore, we also discovered that CS3 has a better prognosis than the patients in CS1 (P = 0.013), CS2 (P = 0.005), and CS5 (P = 0.002), whereas there was no statistical significance compared with the CS4 subgroup. Next, the independent testing for a different distribution of clinical features indicated that these classifications were substantially correlated to overall survival time, age, AJCC stage, Clark level, tumor stage, and metastatic sample type (Table 2).

TABLE 2

Table 2 Summary of clinical features.

To validate the robustness of the cell death clustering, we first applied the “limma” approach to identify the biomarkers for each subtype which are picked from a list of the most up-expressed genes ordered by log₂ fold change (100 biomarkers for each subtype). These biomarkers must satisfy the significance level (adjusted P-value 0.05 and log₂ fold change > 0) and not overlap with any other subtype-identified biomarkers. The heat map visualized the expression of 100 biomarkers for identifying subtypes of melanoma in the TCGA_SKCM cohort (Figure 2D and Supplementary Table S2). Subsequently, based on the template heat map (Figure 2E), nearest template prediction (NTP) was conducted to predict a sample category in the meta-GEO cohort. The Kaplan–Meier survival curve of the predicted five subgroups in the meta-GEO cohort suggested that CS4 was also the best subgroup for prognosis than the others (overall survival, P < 0.001) (Figure 2F). To explore the potential prognostic mechanism in the CS4 subgroup, we mapped the 100 biomarkers in the CS4 subgroup to the STRING website and screened the corresponding proteins. According to the protein–protein interaction (PPI) network analysis (Supplementary Figure S2A), seven hub genes including CD79A, CD79B, CD19, CCR7, CD40LG, SELL, and ZAP70 (degree ≥20) were identified in the CS4 subtype. These hub genes contribute to the activation of primary B-lymphocytes as well as the proliferation of T-cell and cytokine production. The Wilcoxon tests indicated that most of the hub genes (six of seven) were differentially expressed in the TCGA–Genotype Tissue Expression dataset (Supplementary Figure S2B) (24). Furthermore, to prove the results from the database, we also used immunohistochemistry to confirm that these hub genes were also differentially increased in melanoma compared with normal skin samples (Supplementary Figure S2C). Lastly, to test the comparability and repeatability of acquired subtypes across discovery and validation cohorts, the in-group proportion (IGP) statistic will be utilized (25). The values of the IGP scale are from 0 to 1. Most of the IGPs for each subtype in the discovery and validation cohorts were higher than 0.7 (Supplementary Table S3), which indicates a repeatable patient partition for that subtype. Furthermore, several different predictive approaches like partition around medoids classifier, NTP, and consensus ensemble algorithms were performed for accurate classification (21, 26, 27). The Kappa statistics manifested that there is a relatively high consistency among these methods (Figure 3A). Presently, melanomas have traditional classifications such as AJCC stage and Clark level. It is critical to verify the consistency of emerging subtypes with traditional classifications in order to represent the robustness of clustering analysis. Therefore, adjusted mutual information (AMI), Fowlkes–Mallows (FM), Jaccard Index (JI), and Rand Index (RI) were applied to evaluate the similarity between the current subtypes and other well-established subtypes (21). All of these statistics were scaled ranging from 0 to 1, and the higher the number is, the closer the two appraisals are. We observed a moderated agreement in RI, JI, and FM estimation and a slight similarity in AMI calculation. The alluvial diagram also shows the agreement between the AJCC stage, Clark level, and the present subtypes (Figure 3B).

FIGURE 3

Figure 3 Consistency evaluation of different prediction approaches. (A) Consistency heat map using Kappa statistics among different prediction approaches. (B) Agreement of the predicted five subtypes of melanoma with Clark level and tumor stage classification in the cutaneous melanoma of The Cancer Genome Atlas data cohort. (C) Heat map of subtype-specific functional pathways based on upregulated genes for five identified melanoma subtypes.

3.4. Microenvironment, molecular features of subtypes

Firstly, GSEA was conducted for each subtype to identify subtype-specific functional pathways according to their corresponding upregulated genes. Based on the significance threshold (P-value < 0.05 and adjusted P-value < 0.25) and without any overlapped pathways identified for other subtypes, we observed that CS1 mainly enriched in protein synthesis and post-translational modifications such as ribosomal subunit biogenesis, assembly, and methylation. CS2 was significantly involved in energy metabolism and mitochondrial function like oxidative phosphorylation, NADH dehydrogenase complex assembly, and ATP synthesis coupled electron transport. The subtypes of CS3 and CS4 mainly referred to immune response and inflammatory regulation such as positive regulation of receptor binding, interleukin 4 production, and positive regulation of T cell migration. The CS5 subtype was significantly enriched in pathways associated with skin development, keratinization, and cornification (Figure 3C).

Moreover, to explore the relationship between cell deaths and tumor microenvironment, cell death indices (apoptosis, ferroptosis, necroptosis, and pyroptosis) and tumor microenvironment-related predictors (immune score, stromal score, tumor purity, ESTIMATE score, and mRNAsi score) were calculated and plotted. The distribution of cell death indices for each subtype in the TCGA_SKCM cohort is illustrated in Figure 4A. Kruskal–Wallis test determined that melanoma patients in the CS4 subtype had the highest cell death score than the other subtypes, whereas the cell death scores of CS1 and CS2 were located at a much lower level. Similarly, we surprisingly observed that the CS4 subtype had a superior immune score, stromal score, and ESTIMATE score than the other subtypes (Figure 4B), while tumor purity and mRNAsi score in CS4 were the minima among subtypes. The correlation analysis discovered that four kinds of cell death (apoptosis, ferroptosis, necroptosis, and pyroptosis) were not only positively associated with each other but also have a close positive correlation with immune score, stromal score, and ESTIMATE score. There were negative associations between cell deaths with tumor purity and mRNAsi score (Figure 4C). To prove our observation, similar analyses were performed in the meta-GEO cohort. Similar distribution trends and statistical results were observed in cell death indices (Figure 4D) and tumor microenvironment-related predictors (Figure 4E). The circle plot manifested similar correlation coefficients compared with the TCGA_SKCM cohort (Figure 4F).

FIGURE 4

Figure 4 Association between cell deaths and tumor microenvironment. (A) Distribution of four cell death indices (apoptosis, ferroptosis, necroptosis, and pyroptosis) among five identified melanoma subtypes in the cutaneous melanoma of The Cancer Genome Atlas data (TCGA_SKCM) cohort. (B) Distribution of tumor microenvironment-related predictors (immune score, stromal score, tumor purity, ESTIMATE score, and mRNAsi score) among five identified melanoma subtypes in the TCGA_SKCM cohort. (C) Correlation coefficients between cell death indices and tumor microenvironment-related predictors in the TCGA_SKCM cohort. (D) Distribution of four cell death indices among five identified melanoma subtypes in the meta-Gene Expression Omnibus (GEO) cohort. (E) Distribution of tumor microenvironment-related predictors among five identified melanoma subtypes in the meta-GEO cohort. (F) Correlation coefficients between cell death indices and tumor microenvironment-related predictors in the meta-GEO cohort.

We also produced a heat map to comprehensively depict the distinct features of tumor-associated pathways, immune-microenvironment signatures, and expression of immune-checkpoint molecules (Figure 5A). The Kruskal–Wallis test detected that tumor-associated pathways such as cell cycle, m6A regulation, DNA damage response (DDR), mismatch repair, and metabolism hypoxia were positively enriched in the CS1 subtype (Figure 5B). As for the previous report, the immune microenvironment signatures and immune-checkpoint molecules are where the subtypes of CS3 and CS4 had the highest expression values compared to the others (Figure 5B).

FIGURE 5

Figure 5 Landscape of cell deaths for five identified melanoma subtypes. (A) Comprehensive heat map of tumor-associated pathways, immune-microenvironment signatures, and expression of immune-checkpoint molecules among five identified melanoma subtypes. (B) Box plot of tumor-associated pathways, immune-microenvironment signatures, and expression of immune-checkpoint molecules among subtypes. (C) t-SNE plot for melanoma subtype distribution. (D) Graph learning-based dimensionality reduction plot of melanoma subtypes; each color represents a subtype corresponding to the previously defined subtype. (E) The intra-cluster heterogeneity within CS3 subtype, which was further divided into two subgroups according to their location in graph learning. (F) Kaplan–Meier survival curve of subgroups A and B in CS3 subtype. **** means p < 0.0001.

3.5. Pseudotime analysis of subtypes

To better visualize the five clusters based on unsupervised clustering and identify the underlying mechanisms of individual patient distributions, we performed the t-SNE algorithms to investigate and illustrate cluster categorization across samples. The two-dimensional tSNE plot showed that melanoma patients in the TCGA_SKCM cohort can be distinctively partitioned into five clusters (Figure 5C). In addition, the cell death-associated gene expression was subjected to the graph learning-based dimensionality reduction approach. Individual patients were placed on a manifold with sparse tree structures, and the landscape of melanoma for cell deaths was established. Most patients were classified into separate clusters at the end of the branch, such as CS1, CS2, CS4, and CS5 subtypes, which were consistent with our previously described unsupervised clustering (Figure 5D). Individual patient placed in the landscape indicates the general features of cell death in the relevant subtype. In the cell death landscape, for example, the death-active subtype CS4 and death-inactive subtypes CS1 and CS2 were dispersed at the opposite end of the horizontal axis (Figure 5D). As a result, we speculated that the cell death landscape’s vertical axis indicates a different status of cell death for patients.

Moreover, we observed that the samples of the CS3 subtype have scattered distributions on both sides of the vertical axis in the landscape. Therefore, we believed that a significant intra-cluster heterogeneity existed in the CS3 subtype. Based on the vertical axis in the landscape, we successfully divided the CS3 cluster into two subgroups (group A and group B) (Figure 5E). The Kaplan–Meier (KM) survival plots manifested that group A in the CS3 subtype had good survival than group B with log-rank P = 0.046 (Figure 5F). To assess the classification robustness, the CS3 subtype in the meta-GEO cohort and meta-immune response cohort was correspondingly classified into two subgroups. The KM survival curves revealed similar outcomes that group A in the CS3 subtype had a longer survival time than group B regardless of the meta-GEO cohort (Figure 6A) and meta-immune response cohort (Figure 6B). Interestingly, the two subgroups of the CS3 cluster identified by the cell death landscape not only have distinct survival outcomes but also were correlated with different molecular features—for example, patients in group A had a positive enrichment in CD8 T cell effector, antigen processing and presenting machinery (APM), immune-checkpoint, TME score, INFG signature, MHC classes I and II, ICB resistance, and T cell exhaustion (Figure 6C), whereas group B was much more inclined to cancer hallmark pathways like DDR (Figure 6C).

FIGURE 6

Figure 6 Validation of heterogeneity within CS3 subtype and immunotherapy response. (A) Kaplan–Meier survival curve of subgroups A and B in CS3 subtype at the meta-Gene Expression Omnibus (GEO) cohort. (B) Kaplan–Meier survival curve of subgroups A and B in CS3 subtype at the meta-immune response cohort. (C) Box plot of CD8 T cell effector, DNA damage response, antigen processing and presenting machinery, immune-checkpoint, tumor microenvironment score, INFG signature, MHC classes I and II, ICB resistance, and T cell exhaustion between subgroups (A) and (B, D) Immunotherapy response rate of five identified melanoma subtypes in the cutaneous melanoma of The Cancer Genome Atlas data (TCGA_SKCM) cohort. (E) Immunotherapy response rate of five identified melanoma subtypes in the meta-immune response cohort. (F) Rate of immunotherapy response for subgroups A and B in CS3 subtype at the TCGA_SKCM cohort. (G) Rate of immunotherapy response for subgroups A and B in CS3 subtype at the meta-immune response cohort.

3.6. Prediction of immunotherapy and chemotherapy response

Because of the distinct distribution of expression of immune-checkpoint molecules among subtypes, we speculated that our classifications based on cell deaths can be applied to predict immunotherapy response in patients of melanoma. To prove our guesses, the melanoma patients who received immunotherapy in the TCGA_SKCM cohort and the meta-immune response cohort were classified into five categories accordingly. We surprisingly uncovered that the response rate of the CS4 subgroup was relatively higher than the other subgroups no matter if in the TCGA_SKCM cohort or the meta-immune response cohort. Patients in the CS4 subgroup have a 58.33% immunotherapy response in the TCGA_SKCM cohort (Figure 6D). Meanwhile, in the meta-immune response cohort, the CS4 subgroup has the highest response rate (57.14%) than the other subgroups (Figure 6E) Moreover, the CS3 subtype had more diversified cell death and immune characteristics with moderate immunotherapy response. A subdivided analysis indicated that the enrichment score of immune-checkpoints in group A was higher than in group B in the CS3 subtype. Thus, we also explored the immune response rate between group A and group B. We observed that group A had a larger response rate than group B regardless of the TCGA_SKCM cohort and the meta-immune response cohort. In the TCGA_SKCM cohort, patients in group A have a 43.75% response rate versus 29.63% in group B (Figure 6F). At the same time, group A in the meta-immune response cohort has a larger proportion response rate than group B (66.67% vs. 10%) (Figure 6G).

Scientists are struggling to discover new potential compounds for melanoma due to resistance to standard chemotherapeutics (28). To investigate candidate agents that displayed variable effectiveness in cell death gene categorized clusters, the GDSC database was used to test the prediction model (Figure 7A). The AUC values of drug responses were compared within clusters. The GDSC database has 367 compounds in total. Only 219 compounds were evaluated on 34 kinds of melanoma cell lines that were utilized in the analysis. The Kruskal–Wallis tests detected eight compounds, including ACY-1215, CHIR-99021, EHT-1864, ELESCLOMOL, FTI-277, nilotinib, tubastatin A, and TWS-119, that were distinctly affected in cell death gene categorized clusters (Figure 7B). Furthermore, we observed that the AUC values of ACY-1215, EHT-1864, and tubastatin A were significantly lower in cluster 2. Elesclomol and nilotinib were significantly lower in cluster 3. In cluster 5, the AUC values of CHIR-99021 and FTI-277 were much lower than in the other clusters. Furthermore, TWS-119 was much lower in cluster 3. The lower AUC of the drug responses means being more sensitive to drug therapy.

FIGURE 7

Figure 7 Drug identification in cell death-associated gene classified clusters of melanoma. Immunotherapeutic response and potential compounds. (A) Heat map of the differential expression of cell death-associated gene classified clusters in melanoma cells at the Genomics of Drug Sensitivity in Cancer database. (B) Box plot of the area under the curve of ACY-1215, CHIR-99021, EHT-1864, ELESCLOMOL, FTI-277, NILOTINIB, TUBASTATIN A, and TWS-119 among five clusters. (C) Cell viability curves and estimated IC50 values of ACY-1215, tubastatin A, and EHT-1864. (D) Wound healing assay in A375 cells was performed after treatment with ACY-1215 (10.33 μm), tubastatin A (17.77 μm), and EHT-1864 (32.83 μm) at a 48-h recovery period. (E) Flow cytometry analysis of A357 cells which were stained with Annexin V-FITC and propidium iodide after 48 h of ACY-1215 (10.33 μm), tubastatin A (17.77 μm), and EHT-1864 (32.83 μm) treatment.

3.7. Effect of ACY-1215, EHT-1864, and tubastatin A on A357 melanoma cell

We discovered that melanoma patients with death-inactive subtypes are more sensitive to ACY-1215, EHT-1864, and tubastatin A. Therefore, these drugs were selected for cell experiments. Firstly, the CCK8 test was used to measure the cell viability and IC50 values of drugs on A357 melanoma cells. We found that both ACY-1215 and tubastatin A can effectively kill A357 cells with lower concentrations. The IC50 values of ACY-1215, tubastatin, and EHT-1864 were 10.33, 17.7, and 32.83 μm, respectively (Figure 7C). Furthermore, the wound healing assay indicated that ACY-1215, tubastatin, and EHT-1864 can significantly attenuate the cell migration of A357 after incubation for 24 and 48 h (Figure 7D). The effect of ACY-1215, tubastatin, and EHT-1864 on A357 cell apoptosis was studied using flow cytometry. The results show that these drugs can obviously increase A357 cell apoptosis (Figure 7E).

4. Discussion

The prevalence of advanced melanoma has risen steadily over the last few decades. However, treatment regimens have markedly increased in recent years as a result of an improved understanding of the underlying tumorigenesis and signal transduction pathways associated with its pathological process—for example, cancer immunotherapies focused on checkpoint inhibitors have reached significant clinical success. However, its usage is severely limited to the fact that just one-third of patients with most forms of melanoma respond to these inhibitors (29). Several research findings tried to integrate cell death induction with immune-checkpoint inhibitors, resulting in synergistically enhanced antitumor activity, even in checkpoint inhibitor-resistant tumors (17). Under this scenario, researchers proposed a novel notion, immunogenic cell death, which might be induced by radiation, tumor vaccination, immunotherapy, or some form of chemotherapy (30). Immunogenic cell death (ICD) was previously known as immunogenic apoptosis since most kinds of ICD are caused by apoptosis. As the understanding of cell death processes has increased, many nonapoptotic cell deaths have been characterized in recent years. Three well-studied nonapoptotic cell deaths include necroptosis, pyroptosis, and ferroptosis (31–33). As a consequence, in this study, we examined the ICD-associated RNA-seq data of melanoma patients for exploration, independent verification, and comprehensive classification of five reproducible melanoma subtypes in multiple cohort studies, which put a spotlight on cancer classification for frontline therapeutic approaches.

Various gene expression-based biomarkers in the literature can categorize individuals who are at risk and who can benefit from individualized treatments, but their correctness and predictive ability stay restricted (34–38). Nevertheless, the majority of this research relied on the differential study of genomic or transcriptome characteristics rather than biological processes (39). In the current research, we first collected the genes referred to “apoptosis”, “ferroptosis”, “necroptosis”, and “pyroptosis” processes. Moreover, to prove that these genes are consistent with cell death, we conducted GO enrichment analysis and found that these genes are positively activated in cell death-associated signal pathways. Next, we investigated the variations of specific melanoma subtypes based on cell death gene expression and its relationship to tumor genetic mutations as well as the immunological environment. We found that patients in the high-risk subgroups (CS1, CS2, and CS5) with low cell death indices were associated with high Clark levels and tumor stage. These groups also manifested worse overall survival. The substantial connection of cell death with clinical characteristics and survival outcomes suggests that cell death is effective in melanoma patients’ prognosis. Furthermore, we observed that patients in CS3 and CS4 subgroups with good prognoses were significantly enriched in mutations of BRAF and MUC16. Consistently, previous research has proven that BRAF mutation in melanoma patients indicates a higher overall and relative survival rate (40, 41). Several studies indicate that a BRAF mutant could be identified by host immunity and may play a role in anticancer immune responses (42). Erkes et al. recently reported that BRAF inhibitors promoted pyroptosis in anticancer immune responses, pointing to novel treatment approaches for refractory melanoma (43). Furthermore, the mutation of MUC16 was linked to better overall survival in both NSCLC and melanoma (44). Wang et al. have previously proven that MUC16 mutations are mainly enriched in immune-associated pathways and are favorably linked with T-cell activation, which may enhance the prognosis of melanoma patients (45).

Recent research has revealed strong links between cell death and anti-cancer immunity—for example, the tumor-infiltrated CD8+ T cells are thought to induce tumor cell death via perforin-granzyme and Fas-FasL. The new mechanism, however, has illustrated that CD8+ T cells suppress tumor growth by inducing ferroptosis, pyroptosis, and necroptosis, which prompted a reconsideration of the relationship between tumor cell death and immune activation (17). In this study, we observed that apoptosis, necroptosis, pyroptosis, and ferroptosis positively correlated with a highly infiltrated immune microenvironment and highly expressed immune-checkpoint genes, which was consistent with many previous tumor cell death research (46–48). Furthermore, compared with patients in the low-cell-death group, patients in the high-cell-death groups (CS3 and CS4) had a lower expression of tumor-associated molecular pathways and higher enrichment of immune response biological pathways. A variation in the functional immunological network, as previously reported, can cause disruptions in anti-cancer response, immunoediting, and cancer cell escape (49). This heterogeneity of cell death patterns combined with distinct tumor immune microenvironments offers a once-in-a-lifetime chance to develop targeted therapeutics (50). Patients with high infiltration of the immune microenvironment and highly expressed checkpoint target genes are more particularly sensitive to immune-checkpoint inhibitors. We discovered that melanoma patients in CS3 and CS4 subtypes with a high expression of checkpoint target genes are more sensitive to checkpoint inhibitor therapy. However, basic unsupervised clustering analysis for patient subtyping frequently failed to identify intra-cluster interactions and did not reveal the overall structure of the patient distribution (51)—for example, a subtype of CS3 had similar prognosis and clinical outcomes compared with the CS4 subtype. However, the CS3 subtype manifested a more sophisticated heterogeneity of the immune microenvironment and had a lower chromosomal instability than the CS4 subtype. To address these issues, we used graph learning algorithms to identify the tree topologies of cell death profiles in patients, which supplied additional information to clustering analysis and provided fresh insight into the complicated melanoma cell death landscape. We found that the CS3 subtype can successfully cluster into two subgroups (group A and group B) with different survival outcomes, molecular features, and distinct responses of immune-checkpoint inhibitor therapies.

The biological differences between the five categories may signal the necessity for distinct treatment strategies. Based on the GDSC database, we identified eight compounds that were distinctly affected in cell death gene categorized clusters. ACY-1215 and tubastatin A are selective inhibitors of HDAC6, which are required for the proliferation and metastasis of melanoma cells (52, 53). Previous studies have revealed that HDAC6 inhibitors might suppress the growth of a panel of human melanoma cell lines and could be a potential strategy for melanoma therapy, even resolving vemurafenib resistance (53, 54). In our research, we observed that ACY-1215, tubastatin A, and EHT 1864 are more sensitive to the death-inactive subtypes (CS1, CS2, and CS5) of melanoma patients. Cell experiments also confirmed that ACY-1215, tubastatin A, and EHT 1864 can effectively kill A357 cells via apoptosis. Several previous experiments have proven that selective inhibitors of HDAC6 synergistically improve anticancer activity via induction of tumor cell death (55, 56). Elesclomol is an investigational agent that causes oxidative stress, mitochondrial-induced apoptosis in tumor cells, and synergies with taxanes in tumor models (57). In phase II clinical trial of patients with metastatic melanoma, elesclomol with paclitaxel was observed to improve progression-free survival compared with paclitaxel alone (58). The phase II Tasigna Efficacy in Advanced Melanoma trial suggested that nilotinib, a KIT-selective tyrosine kinase inhibitor, may be an effective therapy option for individuals with KIT-mutated advanced melanoma (59). Moreover, CHIR 99021 is a GSK3 kinase inhibitor that can significantly enhance TNF and IFN production, natural cytotoxicity, and antibody-dependent cellular cytotoxicity for effective cancer immunotherapy (60). In conclusion, it is reasonable to speculate that these chemicals may induce various kinds of cell death and can be employed as precise therapies for some death-inactive melanoma subtypes.

5. Conclusion

In summary, we systematically landscaped five cell death subtypes in melanoma that were connected with distinct anti-cancer immunity and found vastly differing patterns in tumor genetic mutation, immune-checkpoint inhibitors, tumor-infiltrating environment, functional orientation, and, most importantly, clinical outcomes. This research gives a conceptual framework for understanding melanoma cell death and immune response and has practical implications for individualized immunotherapy and chemotherapy.

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.

Author contributions

QW was in charge of writing. RW and XW were responsible for the duty of editing. Y-pD was responsible for project administration. All authors contributed to the article and approved the submitted version.

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/fimmu.2022.998454/full#supplementary-material

Supplementary Figure 1 | The Gene Ontology enrichment of cell death- related genes.

Supplementary Figure 2 | Protein–protein interaction (PPI) network of CS4 subtype biomarkers. (A) The PPI network of biomarkers and seven hub genes. (B) Box plot of hub genes between tumor and normal skin groups in the Cancer Genome Atlas–Genotype Tissue Expression dataset. (C) Immunohistochemical images of hub genes in melanoma tissues and normal samples.

Supplementary Table 1 | List of gene sets.

Supplementary Table 2 | Gene biomarkers for identifying the subtypes of melanoma in the cutaneous melanoma of The Cancer Genome Atlas data cohort.

Supplementary Table 3 | The in-group proportion statistic in the cutaneous melanoma of The Cancer Genome Atlas data and meta-Gene Expression Omnibus cohorts.

Abbreviations

TCGA_SKCM, cutaneous melanoma of TCGA data; GEO, Gene Expression Omnibus; ssGSEA, single-sample gene set enrichment analysis; CPI, clustering prediction index; TMB, mutation burden; FGA, fraction genome altered; AUC, area under the curve.

References

1. Siegel RL, Miller KD, Jemal A. Cancer statistics, 2019. CA Cancer J Clin (2019) 69(1):7–34. doi: 10.3322/caac.21551

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Geller AC, Clapp RW, Sober AJ, Gonsalves L, Mueller L, Christiansen CL, et al. Melanoma epidemic: an analysis of six decades of data from the Connecticut tumor registry. J Clin Oncol (2013) 31(33):4172–8. doi: 10.1200/JCO.2012.47.3728

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Robert C, Grob JJ, Stroyakovskiy D, Karaszewska B, Hauschild A, Levchenko E, et al. Five-year outcomes with dabrafenib plus trametinib in metastatic melanoma. N Engl J Med (2019) 381(7):626–36. doi: 10.1056/NEJMoa1904059

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Schachter J, Ribas A, Long GV, Arance A, Grob JJ, Mortier L, et al. Pembrolizumab versus ipilimumab for advanced melanoma: final overall survival results of a multicentre, randomised, open-label phase 3 study (KEYNOTE-006). Lancet (2017) 390(10105):1853–62. doi: 10.1016/S0140-6736(17)31601-X

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Garg AD, Agostinis P. Cell death and immunity in cancer: From danger signals to mimicry of pathogen defense responses. Immunol Rev (2017) 280(1):126–48. doi: 10.1111/imr.12574

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Mattia G, Puglisi R, Ascione B, Malorni W, Care A, Matarrese P. Cell death-based treatments of melanoma:conventional treatments and new therapeutic strategies. Cell Death Dis (2018) 9(2):112. doi: 10.1038/s41419-017-0059-7

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Hu GQ, Yang YJ, Qin XX, Qi S, Zhang J, Yu SX, et al. Salmonella outer protein b suppresses colitis development via protecting cell from necroptosis. Front Cell Infect Microbiol (2019) 9:87. doi: 10.3389/fcimb.2019.00087

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Jorgensen I, Miao EA. Pyroptotic cell death defends against intracellular pathogens. Immunol Rev (2015) 265(1):130–42. doi: 10.1111/imr.12287

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Zhang Z, Shao X, Jiang N, Mou S, Gu L, Li S, et al. Caspase-11-mediated tubular epithelial pyroptosis underlies contrast-induced acute kidney injury. Cell Death Dis (2018) 9(10):983. doi: 10.1038/s41419-018-1023-x

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Dondelinger Y, Delanghe T, Rojas-Rivera D, Priem D, Delvaeye T, Bruggeman I, et al. MK2 phosphorylation of RIPK1 regulates TNF-mediated cell death. Nat Cell Biol (2017) 19(10):1237–47. doi: 10.1038/ncb3608

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Dixon SJ, Stockwell BR. The role of iron and reactive oxygen species in cell death. Nat Chem Biol (2014) 10(1):9–17. doi: 10.1038/nchembio.1416

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Huang S, Cao B, Zhang J, Feng Y, Wang L, Chen X, et al. Induction of ferroptosis in human nasopharyngeal cancer cells by cucurbitacin b: molecular mechanism and therapeutic potential. Cell Death Dis (2021) 12(3):237. doi: 10.1038/s41419-021-03516-y

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Lim B, Greer Y, Lipkowitz S, Takebe N. Novel apoptosis-inducing agents for the treatment of cancer, a new arsenal in the toolbox. Cancers (Basel) (2019) 11(8):1087. doi: 10.3390/cancers11081087

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Lu B, Chen XB, Ying MD, He QJ, Cao J, Yang B. The role of ferroptosis in cancer development and treatment response. Front Pharmacol (2017) 8:992. doi: 10.3389/fphar.2017.00992

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Antonia S, Goldberg SB, Balmanoukian A, Chaft JE, Sanborn RE, Gupta A, et al. Safety and antitumour activity of durvalumab plus tremelimumab in non-small cell lung cancer: a multicentre, phase 1b study. Lancet Oncol (2016) 17(3):299–308. doi: 10.1016/S1470-2045(15)00544-6

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Socinski MA, Jotte RM, Cappuzzo F, Orlandi F, Stroyakovskiy D, Nogami N, et al. Atezolizumab for first-line treatment of metastatic nonsquamous NSCLC. N Engl J Med (2018) 378(24):2288–301. doi: 10.1056/NEJMoa1716948

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Tang R, Xu J, Zhang B, Liu J, Liang C, Hua J, et al. Ferroptosis, necroptosis, and pyroptosis in anticancer immunity. J Hematol Oncol (2020) 13(1):110. doi: 10.1186/s13045-020-00946-7

PubMed Abstract | CrossRef Full Text | Google Scholar

18. 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. doi: 10.1016/j.cell.2016.02.065

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Riaz N, Havel JJ, Makarov V, Desrichard A, Urba WJ, Sims JS, et al. Tumor and microenvironment evolution during immunotherapy with nivolumab. Cell (2017) 171(4):934–949.e916. doi: 10.1016/j.cell.2017.09.028

PubMed Abstract | CrossRef Full Text | Google Scholar

20. 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. doi: 10.1126/science.aad0095

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Lu X, Meng J, Zhou Y, Jiang L, Yan F. MOVICS: An r package for multi-omics integration and visualization in cancer subtyping. Bioinformatics (2020) 14:btaa1018. doi: 10.1093/bioinformatics/btaa1018

CrossRef Full Text | Google Scholar

22. Zeng D, Ye Z, Shen R, Yu G, Wu J, Xiong Y, et al. IOBR: Multi-omics immuno-oncology biological research to decode tumor microenvironment and signatures. Front Immunol (2021) 12:687975. doi: 10.3389/fimmu.2021.687975

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res (2018) 28(11):1747–56. doi: 10.1101/gr.239244.118

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Zeng Ni, Guo C, Wang Y, Li L, Chen Xi, Gao S, et al. Characterization of aging-related genes to predict prognosis and evaluate the tumor immune microenvironment in malignant melanoma. J Oncol (2022) 2022:1271378. doi: 10.1155/2022/1271378

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Kapp AV, Tibshirani R. Are clusters found in one dataset present in another dataset? Biostatistics (2007) 8(1):9–31. doi: 10.1093/biostatistics/kxj029

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Hoshida Y. Nearest template prediction: A single-sample-based flexible class prediction with confidence assessment. PloS One (2010) 5(11):e15543. doi: 10.1371/journal.pone.0015543

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Tibshirani R, Hastie T, Narasimhan B, Chu G. Diagnosis of multiple cancer types by shrunken centroids of gene expression. Proc Natl Acad Sci U.S.A. (2002) 99(10):6567–72. doi: 10.1073/pnas.082099299

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Chasov V, Mirgayazova R, Zmievskaya E, Khadiullina R, Valiullina A, Stephenson Clarke J, et al. Key players in the mutant p53 team: Small molecules, gene editing, immunotherapy. Front Oncol (2020) 10:1460. doi: 10.3389/fonc.2020.01460

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Bagchi S, Yuan R, Engleman EG. Immune checkpoint inhibitors for the treatment of cancer: Clinical impact and mechanisms of response and resistance. Annu Rev Pathol (2021) 16:223–49. doi: 10.1146/annurev-pathol-042020-042741

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Galluzzi L, Vitale I, Warren S, Adjemian S, Agostinis P, Martinez AB, et al. Consensus guidelines for the definition, detection and interpretation of immunogenic cell death. J Immunother Cancer (2020) 8(1):e000337. doi: 10.1136/jitc-2019-000337corr1

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Conos SA, Chen KW, De Nardo D, Hara H, Whitehead L, Nunez G, et al. Active MLKL triggers the NLRP3 inflammasome in a cell-intrinsic manner. Proc Natl Acad Sci U.S.A. (2017) 114(6):E961–9. doi: 10.1073/pnas.1613305114

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Duong BH, Onizawa M, Oses-Prieto JA, Advincula R, Burlingame A, Malynn BA, et al. A20 restricts ubiquitination of pro-interleukin-1beta protein complexes and suppresses NLRP3 inflammasome activity. Immunity (2015) 42(1):55–67. doi: 10.1016/j.immuni.2014.12.031

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Lawlor KE, Khan N, Mildenhall A, Gerlic M, Croker BA, D'Cruz AA, et al. RIPK3 promotes cell death and NLRP3 inflammasome activation in the absence of MLKL. Nat Commun (2015) 6:6282. doi: 10.1038/ncomms7282

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Huang R, Mao M, Lu Y, Yu Q, Liao L. A novel immune-related genes prognosis biomarker for melanoma: associated with tumor microenvironment. Aging (Albany NY) (2020) 12(8):6966–80. doi: 10.18632/aging.103054

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Mei Y, Chen MM, Liang H, Ma L. A four-gene signature predicts survival and anti-CTLA4 immunotherapeutic responses based on immune classification of melanoma. Commun Biol (2021) 4(1):383. doi: 10.1038/s42003-021-01911-x

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Shou Y, Yang L, Yang Y, Zhu X, Li F, Xu J. Determination of hypoxia signature to predict prognosis and the tumor immune microenvironment in melanoma. Mol Omics (2021) 17(2):307–16. doi: 10.1039/D0MO00159G

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Song LB, Luan JC, Zhang QJ, Chen L, Wang HY, Cao XC, et al. The identification and validation of a robust immune-associated gene signature in cutaneous melanoma. J Immunol Res (2021) 2021:6686284. doi: 10.1155/2021/6686284

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Wu Y-L, Hong L-L, Ling Z-N, Hu X-Y, Liu Z, Li P, et al. Golgi scaffold protein PAQR3 as a candidate suppressor of gastric cardia adenocarcinoma via regulating TGF-β/Smad pathway. J Clin Lab Anal (2022) 36(9):e24617. doi: 10.1002/jcla.24617

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Frohlich A, Loick S, Bawden EG, Fietz S, Dietrich J, Diekmann E, et al. Comprehensive analysis of tumor necrosis factor receptor TNFRSF9 (4-1BB) DNA methylation with regard to molecular and clinicopathological features, immune infiltrates, and response prediction to immunotherapy in melanoma. EBioMedicine (2020) 52:102647. doi: 10.1016/j.ebiom.2020.102647

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Heppt MV, Siepmann T, Engel J, Schubert-Fritschle G, Eckel R, Mirlach L, et al. Prognostic significance of BRAF and NRAS mutations in melanoma: a German study from routine care. BMC Cancer (2017) 17(1):536. doi: 10.1186/s12885-017-3529-5

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Wu X, Yan J, Dai J, Ma M, Tang H, Yu J, et al. Mutations in BRAF codons 594 and 596 predict good prognosis in melanoma. Oncol Lett (2017) 14(3):3601–5. doi: 10.3892/ol.2017.6608

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Ilieva KM, Correa I, Josephs DH, Karagiannis P, Egbuniwe IU, Cafferkey MJ, et al. Effects of BRAF mutations and BRAF inhibition on immune responses to melanoma. Mol Cancer Ther (2014) 13(12):2769–83. doi: 10.1158/1535-7163.MCT-14-0290

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Erkes DA, Cai W, Sanchez IM, Purwin TJ, Rogers C, Field CO, et al. Mutant BRAF and MEK inhibitors regulate the tumor immune microenvironment via pyroptosis. Cancer Discovery (2020) 10(2):254–69. doi: 10.1158/2159-8290.CD-19-0672

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Zhang L, Han X, Shi Y. Association of MUC16 mutation with response to immune checkpoint inhibitors in solid tumors. JAMA Netw Open (2020) 3(8):e2013201. doi: 10.1001/jamanetworkopen.2020.13201

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Wang Z, Hou H, Zhang H, Duan X, Li L, Meng L. Effect of MUC16 mutations on tumor mutation burden and its potential prognostic significance for cutaneous melanoma. Am J Transl Res (2022) 14(2):849–62. doi: 10.1001/jamanetworkopen.2020.13201

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Niu X, Chen L, Li Y, Hu Z, He F. Ferroptosis, necroptosis, and pyroptosis in the tumor microenvironment: perspectives for immunotherapy of SCLC. Semin Cancer Biol (2022) 86(Pt 3):273–85. doi: 10.1016/j.semcancer.2022.03.009

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Sheng-Kai H, Chia-Yang L, I-Ling L, Wun-Jyun S, Yih-Fung C, Kai-Chun C, et al. Inflammation-related pyroptosis, a novel programmed cell death pathway, and its crosstalk with immune therapy in cancer treatment. Theranostics (2021) No.18:8813–35. doi: 10.7150/thno.62521

CrossRef Full Text | Google Scholar

48. Zhang Y, Chen Y. Stratification from heterogeneity of the cell-death signal enables prognosis prediction and immune microenvironment characterization in esophageal squamous cell carcinoma. Front Cell Dev Biol (2022) 10:855404. doi: 10.3389/fcell.2022.855404

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Burkholder B, Huang RY, Burgess R, Luo S, Jones VS, Zhang W, et al. Tumor-induced perturbations of cytokines and immune cell networks. Biochim Biophys Acta (2014) 1845(2):182–201. doi: 10.1016/j.bbcan.2014.01.004

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Marzagalli M, Ebelt ND, Manuel ER. Unraveling the crosstalk between melanoma and immune cells in the tumor microenvironment. Semin Cancer Biol (2019) 59:236–50. doi: 10.1016/j.semcancer.2019.08.002

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Li B, Cui Y, Nambiar DK, Sunwoo JB, Li R. The immune subtypes and landscape of squamous cell carcinoma. Clin Cancer Res (2019) 25(12):3528–37. doi: 10.1158/1078-0432.CCR-18-4085

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Liu J, Gu J, Feng Z, Yang Y, Zhu N, Lu W, et al. Both HDAC5 and HDAC6 are required for the proliferation and metastasis of melanoma cells. J Transl Med (2016) 14:7. doi: 10.1186/s12967-015-0753-0

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Woan KV, Lienlaf M, Perez-Villaroel P, Lee C, Cheng F, Knox T, et al. Targeting histone deacetylase 6 mediates a dual anti-melanoma effect: Enhanced antitumor immunity and impaired cell proliferation. Mol Oncol (2015) 9(7):1447–57. doi: 10.1016/j.molonc.2015.04.002

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Peng U, Wang Z, Pei S, Ou Y, Hu P, Liu W, et al. ACY-1215 accelerates vemurafenib induced cell death of BRAF-mutant melanoma cells via induction of ER stress and inhibition of ERK activation. Oncol Rep (2017) 37(2):1270–6. doi: 10.3892/or.2016.5340

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Peng U, Wang Z, Pei S, Ou Y, Hu P, Liu W, et al. ACY-1215 accelerates vemurafenib induced cell death of BRAF-mutant melanoma cells via induction of ER stress and inhibition of ERK activation(Article). Oncol Rep (2017) No.2:1270–6. doi: 10.3892/or.2016.5340

CrossRef Full Text | Google Scholar

56. Li L, Liu F, Huang W, Wang J, Wan Y, Li M, et al. Ricolinostat (ACY-1215) inhibits VEGF expression via PI3K/AKT pathway and promotes apoptosis in osteoarthritic osteoblasts. BioMed Pharmacother (2019) 118:109357. doi: 10.1016/j.biopha.2019.109357

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Albayrak G, Korkmaz FD, Tozcu D, Dogan Turacli I. The outcomes of an impaired powerhouse in KRAS mutant lung adenocarcinoma cells by elesclomol. J Cell Biochem (2019) 120(6):10564–71. doi: 10.1002/jcb.28342

PubMed Abstract | CrossRef Full Text | Google Scholar

58. O'Day S, Gonzalez R, Lawson D, Weber R, Hutchins L, Anderson C, et al. Phase II, randomized, controlled, double-blinded trial of weekly elesclomol plus paclitaxel versus paclitaxel alone for stage IV metastatic melanoma. J Clin Oncol (2009) 27(32):5452–8. doi: 10.1200/JCO.2008.17.1579

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Guo J, Carvajal RD, Dummer R, Hauschild A, Daud A, Bastian BC, et al. Efficacy and safety of nilotinib in patients with KIT-mutated metastatic or inoperable melanoma: Final results from the global, single-arm, phase II TEAM trial. Ann Oncol (2017) 28(6):1380–7. doi: 10.1093/annonc/mdx079

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Cichocki F, Valamehr B, Bjordahl R, Zhang B, Rezner B, Rogers P, et al. GSK3 inhibition drives maturation of NK cells and enhances their antitumor activity. Cancer Res (2017) 77(20):5664–75. doi: 10.1158/0008-5472.CAN-17-0799

PubMed Abstract | CrossRef Full Text | Google Scholar